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bath  of  125  Ar  atoms  at  various  pressures.  The  scaled:  (a)  vibrational  energy,  (b) 
temperature,  and  (c)  pressure.  The  line  coloration:  35  atm,  red;  70  atm,  magenta;  100 
atm,  green;  500  atm,  blue;  900  atm,  yellow;  1300  atm,  dark-green . 112 


Figure  18.  Plots  of  the  decay  of  excess-vibrational  energy  above  the  thermal  limit  (See, 
Fig.  17a)  as  a  function  of  time  on  a  semi-log  scale:  single  exponential  fits  (dashed-black); 
35  atm,  red;  70  atm,  purple;  100  atm,  light-green;  500  atm,  blue;  900  atm,  yellow;  and 
1300  atm  dark-green . 113 

Figure  19.  Same  as  figure  18,  except  scale  is  not  semi-log  and  fitted  curves  are  from  Eq. 
46 . 114 

Figure  20.  Plots  of  the  fitting  amplitudes  and  time  constants  in  Table  X  as  a  function  of 
bath  gas  density.  The  lines  are  intended  to  guide  the  eye:  (a)  the  solid-blue  line 
corresponds  to  the  larger  time  constant  and  the  amplitude  is  matched  with  the  solid-red 
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line.  The  dashed-blue  line  eorresponds  to  the  smaller  time  constant  and  the  amplitude  is 
matched  with  the  dashed-red  line;  (b)  the  results  from  Paul  et  (a)  the  solid-green 
line  corresponds  to  the  larger  time  constant  and  the  amplitude  is  matched  with  the 
solid-yellow  line.  The  dashed-green  line  corresponds  to  the  smaller  time  constant  and  the 
amplitude  is  matched  with  the  dashed-yellow  line.  The  symbols  in  (b)  were  included  to 
emphasize  the  densities  specific  densities  of  the  simulations . 118 

Figure  21  The  low-pressure  time  constants  ki  and  k2  as  a  function  of  pressure  on  a 
smaller  scale  than  that  shown  in  Fig.  19.  The  larger  time  constants  are  plotted  relative  to 
the  axis  on  the  left  and  the  smaller  time  constants  are  plotted  relative  to  the  axis  on  the 
right.  The  solid  and  dashed-black  lines  correspond  to  a  linear  fit  of  the  time  constants 
subject  to  the  constraint  of  zero  energy  transfer  in  the  limit  of  zero  pressure . 120 

Figure  22,  Relaxation  results  as  a  function  of  increasing  numbers  of  Ar  bath  atoms,  (a) 
Scaled  bath  gas  temperature;  250  Ar,  purple  curve;  500  Ar,  green  curve;  750  Ar,  blue 
curve;  1000  Ar,  yellow  curve,  (b)  Same  as  (a),  except  scaled  bath  pressure,  (c) 

Vibrational  energy  decay  on  semi-log  scale  with  the  same  coloration  as  in  (a),  dashed, 
black  lines  correspond  to  fits  to  a  single  exponential  function,  (d)  Vibrational  decay  with 
corresponding  fits  to  Eq.  46 . 122 

Figure  23.  The  variation  of  the  time  constants  and  amplitudes  obtained  in  fitting  Eq.  46 
as  a  function  of  the  number  of  bath  gas  atoms  at  800  K  and  1300  atm.  The  lines  are 
intended  to  guide  the  eye.  The  results  from  the  125  Ar  atoms  at  1300  atm  have  also  been 
included  in  the  plotted  data,  (a)  The  solid-blue  line  correspond  to  the  magnitude  larger 
exponential  and  the  solid-red  line  corresponds  to  the  associated  amplitudes.  The 
dashed-blue  line  corresponds  to  the  magnitude  of  smaller  exponential  and  the  dashed-red 
lines  correspond  to  the  associated  amplitudes,  (b)  Plots  of  the  magnitude  of  the  time 
constants  on  an  expanded  scale;  larger  time  constants  are  in  blue  and  are  relative  to  the 
left  y-axis  and  the  smaller  time  constants  are  in  yellow  and  are  relative  to  the  right  axis 
. 124 
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ABSTRACT 

The  first  study  of  this  dissertation  is  focused  on  studies  of  the  classical  dynamics 
and  associated  rates  of  isomerization  and  dissociation  of  isolated  HO2.  Comparisons  are 
made  between  the  three  potential  energy  surfaces  (PESs)  used  in  these  studies.  The 
intramolecular  vibrational  energy  redistribution  (IVR)  at  energies  above  and  below  the 
threshold  of  isomerization  is  slow,  especially  for  0-0  stretch  excitations,  consistent  with 
the  regularity  in  the  surfaces-of-section.  The  slow  IVR  rates  lead  to  mode-specific  effects 
that  are  prominent  for  isomerization  and  modest  for  unimolecular  dissociation  to  H  +  O2. 
Even  with  statistical  distributions  of  initial  energy,  slow  IVR  rates  result  in 
bi-exponential  decay  for  isomerization,  with  the  slower  rate  correlated  with  slow  IVR 
rates  for  0-0  vibrational  excitation.  The  calculated  IVR  results  for  all  three  PESs  are 
reasonably  well  represented  by  an  analytic,  coupled  three-mode  energy  transfer  model. 

The  second  study  of  this  thesis  is  focused  on  the  effects  of  pressure  on  the 
relaxation  of  the  HO2  embedded  in  a  dense  gas  environment.  A  method  of  simulating  the 
radical  in  an  argon  bath  is  proposed  and  validated.  The  time-dependent  decay  of 
vibrational  energy  is  found  to  be  bi-exponential  for  all  of  the  simulated  pressures.  The 
relaxation  rates  at  low  pressures  extrapolate  poorly  to  the  high-pressure  results,  with  a 
turnover  in  the  rate  constants  occurring  at  intermediate  pressures.  The  effects  of  finite 
size  on  the  simulation  are  investigated.  Comparisons  to  studies  with  similar  findings  and 
additional  considerations  for  understanding  this  behavior  are  discussed. 
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I.  Introduction 

The  interplay  of  theory  and  experimentation  has  been  a  erueial  engine  advaneing 
the  understanding  the  mechanisms  of  unimolecular  reactions/  The  progression  of 
theories  describing  unimolecular  reactions  are  discussed  in  review  articles  and  in  several 
monographs/  The  Lindemann  mechanism"^  is  the  foundation  on  which  more  recent 
theories  are  built  upon.  An  example  of  this  mechanism  is  the  following 


A  +  M  :^A*+M 

(1) 

A^^B, 

(2) 

where  ^4  is  a  reactant,  ^4  *  is  a  vibrationally  excited  reactant,  5  is  a  product  and  M  is  a 
buffer  gas.  The  gist  of  the  theory  is  that  a  chemical  species  gains  and  loses  energy  by 
collisions  with  the  surrounding  medium  which  can  lead  to  a  reaction  occurring.  At  low 
pressures  this  is  a  bimolecular  process,  and  at  higher  pressures  the  process  in  Eq.  2 
becomes  the  rate  determining  step  of  the  mechanism.  However,  this  mechanism  only 

-5 

qualitatively  predicts  the  correct  behavior  of  the  reaction  rate  as  a  function  of  pressure. 
The  issue  is  that  the  internal  state  of  the  reactant  must  be  addressed  for  better  accuracy. 
The  theories  of  Hinshelwood,  followed  by  those  of  Ramsberger,  Rice,  and  Kassel 

o 

(RRK)  and  Slater  all  attempted  to  resolve  the  issue  of  how  the  internal  energy  of  the 
chemical  species  affects  the  rate  of  reaction.  These  earlier  theories  provided  a  foundation 
for  the  extensions  of  RRK  proposed  by  Marcus,^  which  has  become  known  as  RRKM 
theory.  RRKM  theory  extends  the  ideas  of  RRK  by  incorporating  a  quantum-statistical 
mechanical  description  of  the  energization  of  a  reactant  species  and  transition  state 
theory  ’  to  describe  the  conversion  of  the  energized  reactant  to  a  product  state.  A  key 
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assumption  is  that  localization  of  energy  in  a  vibrational  degree  of  freedom  is 
unimportant  as  long  as  the  intramoleeular  redistribution  of  energy  (IVR)  amongst  the 
vibrational  modes  is  fast  relative  to  the  timescale  of  the  reaction.  An  expectation  of  the 
assumption  of  fast  IVR  is  that  unimolecular  reactions  should  occur  in  a  probabilistic 
manner  and  where  the  reaetion  of  a  ehemical  species  is  independent  of  its  initial  state.  A 
brief  aside  is  neeessary  to  discuss  the  nature  of  energy  transfer  amongst  vibrational 
modes. 


The  first  eomputational  studies  of  nonlinear  dynamics  performed  by  Fermi,  Pasta, 
and  Ulam^°  (FPU)  on  the  equipartion  of  energy  in  a  one-dimensional  anharmonie  chain 
of  oscillators  yielded  results  that  surprised  the  authors.  The  simulation  consisted  of 
exciting  a  single  vibrational  mode  in  a  one-dimensional  chain  of  oscillators  and 
propagating  the  equations  of  motion.  The  intent  was  to  study  how  nonlinear  forces 
affeeted  the  transfer  of  energy  to  the  higher  vibrational  modes  of  the  system.  However, 
they  found  that  on  a  long  time-seale,  there  was  non-ergodic  behavior  with  an  almost 
eomplete  reloealization  of  energy  in  the  excited  mode,  which  contradicted  the  expeeted 
equilibration  of  energy  amongst  the  modes  of  the  system.  Additional  work  by  Tuck  and 
Menzel^^  showed  that  this  long-time  periodieity  in  the  FPU  work  was  in  fact  real  and  was 
not  due  to  numerieal  error.  Ford  showed  the  lack  of  relaxation  in  the  FPU  problem  was 
due  to  the  absence  of  resonanees.  The  effeet  of  resonance  conditions  on  ergodic  behavior 
in  trajeetories  is  generalized  in  the  Kolmogorov,  Amol’d,  and  Moser  theorem, which 
states  that  under  a  small  perturbation  periodic,  or  quasi-periodic  motion  in  phase  space 
will  persist  unless  resonant  eonditions  oceur.  This  aside  on  nonlinear  dynamies  ties 
direetly  to  a  question  that  has  been  posed  by  others:  Does  the  flow  of  energy  amongst  the 
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vibrational  modes  of  a  molecular  system  affect  the  rate  of  reaction?^*^’’^’'^  The  possible 
connection  between  non-statistical  dynamics  and  the  unimolecular  reactivity  motivated 
early  classical  trajectory  studies  of  model  tri-atomics  by  Bunker/"^  He  found  that 
frequency  mismatches  between  vibrational  modes  due  to  mass  differences  could  cause 
non-random  lifetimes  in  the  dissociating  species.  Bunker  and  Hase'^  further  defined  non- 
statistical  lifetimes  by  differentiating  between  effects  due  to  localization  of  energy  and 
effects  due  to  the  phase  space  of  chemical  species  being  metrically  decomposable. 

Some  expectations  of  the  dynamics  of  HO2  can  be  guided  by  observations  that 
others  have  made  while  studying  H2O2.  There  are  numerous  studies^^  which  have 
furthered  our  understanding  of  this  system;  however,  for  brevity  we  focus  on  a  few  of  the 
key  conclusions.  First,  the  motion  of  the  HO  is  relatively  decoupled  from  the  other 
degrees  of  freedom  of  the  system.  Second,  there  is  evidence  of  slow  energy  flow  out  of 
an  excited  HO  mode  that  has  been  prepared  in  an  overtone  state.  This  slow  energy  flow 
could  contribute  to  non-exponential  population  decay  of  a  sample.  Finally,  participation 
of  the  lower-frequency  modes  and  rotational  motion  along  an  axis  defined  by  the  0-0 
bond  has  the  effect  of  increasing  the  rate  of  IVR.  While  the  addition  of  a  hydrogen  atom 
to  HO2  leads  to  considering  the  motion  of  three  more  internal  degrees  of  freedom,  the 
vibrational  frequencies  of  these  two  species  share  similarities,  so  one  would  expect  some 
parallels  between  the  intramolecular  dynamics  of  these  two  species.  More  specifically, 
slow  energy  transfer  between  the  high  frequency  and  low  frequency  modes. 

The  first  part  of  this  dissertation  is  tied  to  the  simplified  reaction  shown  in  Eq.  2, 
namely  how  does  the  unimolecular  reaction  of  HO2  occur?  Are  the  dynamics  in 
agreement  with  statistical  expectations,  with  fast  energy  transfer  among  the  vibrational 
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modes?  What  role  does  isomerization  play  in  the  energy  transfer?  Are  there  signatures 
of  mode  selectivity  in  the  unimolecular  dissociation  of  the  radical?  These  questions  are 
the  subject  of  the  research  in  Chapter  III  of  this  dissertation. 

The  second  part  of  this  dissertation  is  devoted  to  the  study  of  the  collisional 

deactivation  of  the  vibrationally  excited  species  shown  in  Eq.  1.  The  isolated  binary 

collision  approximation  (IBC)  is  a  common  assumption  made  in  modeling  the 

deactivating  effect  of  collisions  on  a  vibrationally  excited  species.  Under  this 

approximation,  the  time-dependent  deactivation  process  is  treated  as  a  series  of 

uncorrelated  collisions  with  the  bath  gas.  This  simplified  model  allows  for  the 

calculation  of  the  relaxation  time  of  a  solute  from  the  product  of  an  energy  transference 

probability  with  a  solute-bath  gas  collision  frequency.  The  energy  transference 

probability  is  typically  based  on  gas-phase  results,  and  the  collision  frequency  is  taken 

from  a  hard-sphere  model.  A  general  result  of  the  IBC  approximation  is  a  linear 

dependence  of  the  solute  relaxation  time  constants  with  the  bath  density.  Evidence  for 

1 8 

the  general  applicability  of  IBC  is  mixed,  with  favorable  agreement  for  some  systems, 
whereas  anomalous  behavior  is  seen  in  others.  Given  the  importance  of  HO2  in 
pressure-dependent  combustion,  we  investigate  the  relaxation  of  a  vibrationally  excited 
radical  in  a  dense  Ar  bath.  We  also  examine  the  applicability  of  the  IBC  approximation 
under  conditions  of  extreme  pressure.  These  investigations  are  the  focus  of  chapter  IV. 

The  development  of  the  code  used  in  simulating  HO2  led  to  preliminary  work  on 
general  methods  and  considerations  of  molecular  dynamics  (MD).  This  work  does  not 
directly  fit  in  other  chapters,  but  it  was  useful  in  motivating  the  work  on  HO2  and  is 
discussed  in  chapter  II.  The  conclusions  of  the  dissertation  are  given  in  chapter  V. 
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II.  Considerations  in  Trajectory  Calculations 

While  the  methods  in  an  MD  simulation  ean  beeome  quite  involved,  the  teehnique 
can  be  broken  into  a  simple  set  of  steps:  (1)  definition  of  the  potential  energy  surface;  (2) 
selection  of  initial  conditions;  (3)  integrating  the  equations  of  motion;  (4)  sampling  and 
analysis  of  physical  properties  and  end-test  criteria.  These  steps  are  shown  below  for  the 
simple  case  of  trajectories  integrated  on  the  Henon-Heiles  potential.  Extension  of  these 
steps  to  more  complicated  simulations  is  also  discussed.  The  discussion  of  this  chapter 
sets  the  stage  for  the  research  conducted  in  completing  this  dissertation. 

A.  The  Henon-Heiles  Hamiltonian:  An  Illustrative  Example 

Defining  a  potential  energy  surface  (PES)  is  the  first  step  in  an  MD  simulation. 

The  purpose  of  the  PES  is  to  represent  the  interactions  of  the  nuclei  that  compose  the 

chemical  species  being  studied.  Originally  the  Henon-Heiles  potential  was  developed  as 

a  model  for  understanding  cosmological  motion,  but  the  simple  form  and  complex 

2 1 

dynamics  of  this  potential  has  resulted  in  its  use  in  understanding  chemical  dynamics. 
The  Henon-Heiles  potential  is  the  following: 

=  +y^)+x^y-^y\  (3) 

where  x  and  y  are  Cartesian  coordinates.  This  potential  is  two  harmonic  oscillators  with 
cubic  anharmonic  coupling  and  anharmonicity  in  the  y  coordinate.  An  illustration  of  the 
potential  can  be  seen  in  Eig.  E  The  minimum  of  the  potential  is  the  central  blue  region. 
The  coloration  of  the  figure  changes  with  progressively  higher  potential  energy  from 
green  to  yellow  and  finally  red.  The  extremely  high-energy  regions,  which  are  of  little 
dynamical  importance,  are  shown  as  white.  The  potential  has  three-fold  symmetry  with 
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the  exit  channels  diverging  to  infinity  both  spatially  and  energetically.  The  threshold  for 
a  dissociative  trajectory  is  ~l/6  energy  units  and  trajectories  with  energy  beneath  that 
energy  threshold  display  behavior  characteristics  that  range  from  quasi-periodic  to 
irregular. 


Figure  1  Henon-Heiles  potential  plotted  in  Cartesian  coordinates.  The  black  lines 
indicate  planes  of  intersection  that  were  used  in  calculating  surfaces-of-section. 


The  initial  conditions  for  these  illustrative  tests  were  selected  with  a  simple  Monte 
Carlo  scheme.  The  first  step  involves  selecting  two  pseudo-random  numbers  uniformly 
distributed  on  [-1,1]  as  the  initial-Cartesian  coordinates.  The  potential  is  evaluated  at 
these  points  to  determine  if  they  lie  below  the  specified  total  energy.  If  the  potential 
energy  is  above  the  specified  threshold,  the  points  are  then  discarded  and  a  new  set  of 
points  is  sampled;  this  process  is  repeated  until  the  energy  criterion  is  satisfied.  Once  the 
energy  criterion  is  met,  then  a  set  of  conjugate  momenta  is  selected  that  satisfies  T  =  H-V 
where  T  is  the  kinetic  energy,  H  is  the  total  energy  and  V  is  the  potential.  The  kinetic 
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energy  is  then  divided  between  the  two  oscillators  using  and 

Py  -  yl2T(\-(^)  ,  where  C  is  a  pseudo-random  number  uniformly  distributed  on  [0,1],  Px 

is  the  momentum  in  the  x  direction  and  Py  is  the  momentum  in  the  y  direction.  The  signs 
of  the  momenta  are  randomized  and  a  trajectory  is  initiated.  This  trajectory  is  then 
integrated  for  some  desired  length  of  time,  or  until  a  termination  criterion,  such  as 
dissociation,  is  satisfied.  The  trajectory  terminates  when  a  predefined  time  or  event 
occurs  and  the  initial  conditions  of  the  next  trajectory  are  selected.  The  process  is 
repeated  until  a  group,  or  ensemble,  of  trajectories  is  completed.  During  the  simulation, 
the  trajectory  can  be  analyzed  on  the  fly  with  ensemble  averages  being  calculated  at  the 
end  of  the  simulation,  or  phase  points  can  be  stored  for  subsequent  post-simulation 
analysis.  Collecting  the  full  phase  space  information  has  the  advantage  that  new  analyses 
can  be  considered  without  rerunning  trajectories.  The  main  disadvantage  of  storing  the 
entire  phase  space  information  of  a  simulation  is  the  dependence  on  having  sufficient 
disk  space  available  and  preventing  data  corruption.  The  on-the-fiy  analysis  has  the 
advantage  of  requiring  minimal  disk  space  and  data  corruption  issues,  but  new  analysis 
can  require  rerunning  simulations.  The  complex  dynamics  of  the  Henon-Heiles  potential 
has  resulted  in  its  extensive  study,  particularly  with  regard  to  the  transition  from  quasi- 
periodic  dynamics  to  irregular  dynamics  where  the  phase  space  of  the  system  is  explored 
in  a  much  more  chaotic  manner.  Figure  2  contains  two  representative  trajectories  which 
are  representative  of  quasi-periodic  (frame  a)  and  irregular  dynamics  (frame  b). 
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Figure  2.  Representative  trajeetories  of  quasi-periodic  and  irregular  trajeetories:  (a) 
Quasi-periodie  trajectory  with  a  total  energy  of  0.05  energy  units;  (b)  irregular  trajectory 
with  a  total  energy  of  0.16  energy  units;  (c)  and  (d)  the  X  surface-of-section;  (e)  and  (f) 
the  Y  surface-of-section. 

Now  imagine  attempting  to  classify  an  ensemble  of  the  trajectories.  Examining 
plots  of  individual  trajectories  would  be  inefficient,  and  understanding  the  behavior  of  the 
ensemble  would  be  difficult.  A  slight  change  in  perspective  can  provide  some  clarity 
shown,  as  in  Fig.  2c,  2d,  2e  and  2f  by  using  Poincare  surfaces-of-section  (SOS).  An  SOS 
plot  is  calculated  by  first  selecting  a  plane  defined  by  one  of  the  coordinates  in 
configuration  space  and  recording  the  second  coordinate  and  conjugate  momentum  when 
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a  trajectory  crosses  the  surface  in  one  direction.  A  surface-of-section  is  a  common 
mapping  technique  in  the  analysis  of  dynamical  systems  with  low  dimensionality. 
Higher  dimensional  systems  often  require  more  involved  analysis.  An  SOS  plot 
provides  a  striking  means  of  showing  recurrence,  or  periodicity,  in  an  group  of 
overlapping,  complicated  trajectories.  Figure  3  shows  how  the  surfaces-of-section  for  an 
ensemble  of  trajectories  change  as  the  energy  of  the  system  increases.  As  the  energy  of 
the  system  increases,  the  amount  of  accessible  phase  space  volume  grows  and  the  regions 
of  structure  begin  to  disappear  as  the  trajectories  become  more  irregular. 
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Figure  3.  Surfaces-of-section  from  dynamics  on  the  Henon-Heiles  potential  that  show  the 
progression  of  phase  spaee  expansion  and  exploration  with  increasing  energy.  The 
trajectories  were  integrated  for  1000  time  units.  Eaeh  frame  is  the  result  of  intersections 
calculated  from  ensembles  of  100  trajeetories  per  frame.  The  planes  of  interseetion  are 
defined  in  Fig.  1;  (a)  and  (b)  the  total  energy  is  0.042;  (c)  and  (d)  the  total  energy  is 
0.084;  (e)  and  (f)  the  total  energy  is  0.125;  (g)  and  (h)  the  total  energy  is  0.166. 
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B.  Potential  Energy  Surfaces 

While  a  full  quantum  mechanical  solution  for  the  motion  of  both  the  electrons  and 
the  nuclei  is  desired,  it  is  often  cost  prohibitive.  But,  to  a  good  approximation,  an 
acceptable  solution  can  often  be  obtained.  The  good  approximation  is  otherwise  known 
as  the  Born-Oppenheimer  approximation,  where  the  general  idea  is  that  the  motion  of 
the  electrons  can  be  separated  from  the  motion  of  the  nuclei  by  virtue  of  the  vastly 
differing  masses  of  these  two  subatomic  particles.  In  other  words,  the  motion  of  the 
electrons  is  so  fast  that  one  can  conceive  the  electrons  instantaneously  adjusting  for  the 
motion  of  the  nuclei,  and  that  the  nuclei  move  on  a  surface  that  is  defined  by  energies 
obtained  by  solving  the  electronic  structure  of  the  system.  However,  given  the 
computational  cost  of  electronic  structure  calculations,  it  is  common  practice  to  construct 
model  systems  with  analytic  functions  describing  the  interactions  between  the  nuclei. 
These  analytic  functions  can  be  parameterized  to  reflect  the  known  information  on  the 
chemical  system,  such  as  the  thermochemistry  and  vibrational  frequencies  of  the  system. 

The  analytic  functions  parameterized  with  experimental  and  ab  initio  data  may 
not  be  sufficient  to  represent  the  complex  topography  of  the  true  PES.  But  there  are  two 
different  means  of  obtaining  a  more  realistic  PES.  With  the  improvement  in  computing 
power  and  electronic  structure  codes,  there  has  been  development  of  a  method  known  as 
“Direct  Dynamics.”  With  Direct  Dynamics,  the  electronic  structure  problem  is  solved 
on  the  fly  during  a  trajectory.  The  primary  issue  with  this  technique  is  the  sheer 
computational  cost  of  solving  the  electronic  structure  problem  for  tens  of  thousands  of 
force  evaluations  in  a  single  trajectory.  This  cost  is  compounded  by  the  need  to  simulate 
hundreds,  if  not  thousands,  of  trajectories  to  obtain  converged  averages.  The  alternative 
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to  direct  dynamics  is  to  save  the  output  of  electronic  structure  calculations  and  interpolate 

26 

it  using  a  flexible  function  or  polynomial.  However,  interpolation  of  low-symmetry 
systems  that  contain  four  or  more  atoms  suffers  from  the  “curse  of  dimensionality,”  in 
that  a  large  number  of  data  points  will  be  required  to  fit  a  high-dimension  PES  accurately. 
Another  issue  with  interpolation  is  that  the  computational  cost  grows  with  the  complexity 

98 

of  the  function  being  used  to  fit  the  surface. 

In  the  context  of  this  dissertation  and  the  HO2  radical,  we  are  fortunate  to  have 
access  to  interpolated  potentials  that  were  developed  by  others. We  use  these  PESs  to 
study  the  phenomena  in  which  we  are  interested,  which  are  the  dynamics  and  reactivity 
of  the  HO2  radical.  The  description  of  the  potential  energy  surfaces  of  this  radical  used  in 
this  dissertation  is  provided  in  Ch.  III. 

C.  Selection  of  Initial  Conditions 

The  selection  of  initial  conditions  for  trajectories  is  very  much  dependent  on  the 
statistical  ensemble  and  what  one  wishes  to  learn  from  the  trajectories.  Whatever  the 
phenomenon  of  interest,  the  Monte  Carlo  method  ’  ’  is  a  common  means  for  selecting 
initial  conditions  Typically,  one  uses  a  pseudo-random  number  generator  to  select 
random  numbers  distributed  on  an  interval  between  0  and  1 .  These  numbers  can  then  be 
weighted,  or  scaled,  to  the  variables  needed  to  describe  the  initial  state  of  the  system.  Eor 
NPT  (constant  temperature,  pressure,  and  number  of  particles),  and  NVT  (constant 
temperature,  volume,  and  number  of  particles)  simulations,  the  initial  conditions 
generally  consist  of  selecting  initial  coordinates  and  momenta  and  then  running  a 
sufficiently  long  trajectory  with  a  thermostat  or  barostat  to  obtain  the  desired  physical 
state.  The  sort  of  initial  conditions  that  are  of  greatest  use  in  studying  unimolecular 
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dynamics  are  ones  that  eorrespond  to  the  NVE  ensemble  (eonstant  energy,  volume  and 
number  of  partieles).  This  ensemble  is  attraetive  because  there  are  general  eonservation 
prineiples:  energy,  total  angular  momentum,  and  linear  momentum,  whieh  ean  be  used  to 
assess  the  aceuraey  of  the  simulation.  But  most  important,  the  equations  of  motion  are 
not  modified  as  they  would  be  with  a  thermostat  or  a  barostat,  so  in  the  elassieal  limit 
where  quantum  effeets  are  assumed  to  be  minimal,  a  trajeetory  ealeulated  on  an  aecurate 
PES  should  reproduee  the  dynamies,  and  the  ensemble  averages  ealeulated  by  trajeetories 
would  agree  with  the  expeetation  values  derived  from  quantum  mechanies.  There  has 
been  eonsiderable  investment  of  time  and  effort  in  the  development  of  the  algorithms 
used  to  seleet  initial  eonditions  that  eorrespond  to  the  NVE  ensemble,  so  we  will  foeus  on 
the  ones  most  relevant  to  this  study. 

1.  Mode-Specific  Initial  Conditions 

The  ability  to  loealize  energy  in  a  partieular  degree  of  freedom  is  of  interest 
inasmuch  as  it  provides  a  means  of  studying  how  energy  relaxes  amongst  the  various 
degrees  of  freedom.  These  methods  have  been  diseussed  in  great  detail  in  numerous 
reviews  on  elassieal  simulations  and  are  only  briefly  deseribed  here.  The  simplest  form 
of  approximating  the  vibrational  degrees  of  freedom  is  to  picture  them  as  normal  modes. 
The  normal  mode  approximation  provides  a  general  means  of  separating  the  vibrational 
motion  into  diserete  parts.  However,  this  method  does  have  limitations:  the  vibrational 
motion  is  assumed  to  be  harmonic  and  there  is  no  aeeounting  for  the  anharmonieity  of  the 
PES.  The  form  of  a  normal  mode  is 

(t)  =  A-  cos(/l'^^t  +  ,  (4) 
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where  qt  is  the  normal  mode,  At  is  the  amplitude,  A  is  the  foree  eonstant  and  ^ is  the 
phase  of  the  oseillator.  The  first  step  in  setting  the  mode-speeifie  initial  eonditions  is  to 
solve  the  eigenvalue  equation  of  the  normal  modes  for  the  eigenvalues  and  veetors.  This 
is  done  by  ealeulating  the  Hessian,  or  seeond  derivative  matrix  of  the  potential  at  a 
stationary  point  on  the  PES,  sueh  as  a  minimum  on  the  surfaee.  Detailed  discussions  of 

-5  1 

normal  mode  analysis  and  spectroscopy  can  be  found  in  a  classic  text.  The  eigenvalue 
equation  is  then  solved,  typically  with  readily  available  codes.  The  eigenvector  matrix 
is  used  in  the  linear  transformation  between  normal  coordinates  and  Cartesian 
coordinates  during  the  selection  of  initial  conditions,  but  can  also  find  use  in  the  analysis 
of  the  trajectory.  The  selection  of  initial  conditions  begins  by  calculating  the  total  energy 
summed  from  the  individual  mode  contributions.  Then  the  phase  of  each  oscillator  is 
sampled  independently  using 

^  =  +  (5) 

where  qe  is  the  equilibrium  value  of  the  normal  mode  and  qt  is  the  inner  turning  point 

and  is  a  pseudo-random  number.  The  potential  energy  of  the  surface  is  then  calculated 
with  the  sampled  phase  and  is  accepted  if  the  potential  energy  of  the  displacement  is  less 
than  the  specified  energy  of  the  mode.  If  the  potential  energy  is  greater  than  the  specified 
mode  energy  then  displacement  is  rejected  and  the  phase  is  resample  until  an  acceptable 
displacement  is  found.  After  independently  sampling  the  phase  of  all  of  the  modes,  the 
normal  mode  velocities  and  displacements  are  added  to  the  stationary  point  at  which  the 
normal  mode  analysis  was  done.  The  resultant  sum  of  the  coordinate  and  momentum 
vectors  may  result  in  a  non-zero  angular  momentum,  which  Hase  et  al.  showed  can  be 
iteratively  removed  to  enforce  a  requirement  of  zero  angular  momentum  if  it  is  desired. 
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The  use  of  normal  mode  exeitations  has  found  eonsiderable  use  in  studying  IVR.  A 
eomplementary  method  of  loealizing  energy  in  a  partieular  mode,  similar  to  normal  mode 
exeitation  ean  be  done  using  a  loeal  mode  model, where  the  exeitation  energy  is 
projeeted  along  a  bond  veetor  with  the  assumption  of  the  potential  interaetion  being 
harmonie,  or  Morse-like.  The  advantage  of  a  loeal  mode  exeitation  is  that  the 
anharmonieity  of  the  potential  in  the  loeal  mode  is  taken  into  aeeount  in  the  selection  of 
the  initial  conditions. 

2.  Statistical  Initial  Conditions 

The  limitation  of  assuming  a  normal  or  local  mode  picture  has  led  to  development 
of  alternative  means  of  selecting  initial  conditions  for  the  microcanonical  ensemble.  The 
question  of  efficient  and  accurate  selection  of  initial  conditions  for  the  microcanonical 
ensemble  has  been  extensively  studied  by  others  The  initial  conditions  of  the  HO2 
radical  were  selected  using  efficient  microcanonical  sampling  (EMS)  with  no  restriction 
to  angular  momentum  or  with  the  angular  momentum  restricted  to  zero.  The  procedures 
consists  of  a  random  walk  where  the  Cartesian  coordinates  of  the  radical  are  randomly 
changed 

tjnew  (6) 

where  qnew,  and  qoid  are  the  “old”  and  “new”  x,  y,  or  z  component  of  the  Cartesian 
coordinates  of  the  atoms  and  C  is  a  pseudo-random  number  on  the  range  0  to  1  and  zl  is  a 
maximum  displacement.  The  acceptance  of  a  new  displacement  with  no  restriction  to 
angular  momenta  is  assigned  the  following^"^*^’’^  weight 

r(£,J)  =  [£-F(?)]'“-=>'3  (7) 
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and  the  corresponding  weight  for  J  =  0 

W{E,J)  =  ,  (8) 

where  la,  h,  4  are  the  principal  moments  of  inertia,  E  is  the  desired  energy,  V  is  the 
current  potential  energy,  and  N  is  the  number  of  atoms  in  the  system.  A  trial  move  is 
accepted  if  the  weight  of  the  new  state  is  greater  than  weight  of  the  current  configuration. 
If  the  trial  weight  is  less  than  or  equal  to  the  current  weight,  then  a  ratio  of  the  trial 
weight  and  the  current  weight  is  compared  to  a  pseudo-random  number  sampled  from  a 
uniform  distribution  of  [0,1].  If  the  ratio  of  weights  is  greater  than  the  selected  pseudo¬ 
random  number  the  trial  move  is  then  accepted;  otherwise,  it  is  rejected  and  a  new  trial 
displacement  is  sampled.  At  the  end  of  the  random  walk,  a  set  of  velocities  is  selected 
and  scaled  to  ensure  that  the  desired  initial  internal  energy  Eo  is  obtained.  The 
momentum  scaling  is 


where  El  is  the  Hamiltonian,  and  Pq  and  are,  respectively,  the  old  and  new  x,  y,  or  z 
momenta  of  the  atoms  in  the  radical.  The  final  geometric  configuration  of  system  at  the 
end  of  the  Markov  walk  was  used  as  the  starting  point  for  subsequent  random  walks  in 
the  ensemble. 

D.  Integration  of  the  equations  of  motion 

After  selecting  the  initial  conditions,  the  next  task  is  integration  of  the  equations 
of  motion.  Typically  a  closed  form  solution  to  the  equations  of  motion  does  not  exist  and 
a  numerical  solution  is  required.  Numerical  integration  is  a  rich  and  varied  field  and  has 
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been  deeply  studied  in  the  eontext  of  MD  simulations.  ’  While  no  development  of 
integration  methods  oeeurred  in  this  dissertation,  eoding  and  testing  of  existing 
algorithms  did  happen,  so  it  is  worthwhile  to  eonsider  the  range  of  integration  algorithms 
available  and  the  respeetive  strengths  and  weakness  of  these  algorithms. 

In  MD  simulations,  the  equations  of  motion  solved  are  typieally  Hamilton’s 
equations  sinee  these  equations  are  valid  in  any  eoordinate  frame  and  the  solving  a  set  of 
first-order  equations  is  generally  less  eomplieated  than  solving  a  set  of  seeond-order 
equations.  Hamilton’s  equations  are  the  following 


and 


P 


dH 

dq 


(10) 


dH 

dp 


(11) 


where  p  and  q  are  the  time  derivatives  of  the  momenta  and  eoordinates,  respeetively, 
and  H  is  the  Hamiltonian,  which  is  a  function  of  p  and  q.  Typically,  the  integration 
coordinates  are  selected  to  be  in  the  fixed  lab-frame  Cartesian  coordinates,  a  choice  that 
simplifies  the  expression  of  the  kinetic  energy.  This  expression  can  be  complicated  if  the 
system  is  described  in  terms  of  in  generalized  coordinates.  The  issue  of  the  potential 
energy  being  a  function  of  internal  coordinates  is  dealt  with  by  using  the  chain  rule  to 
transform  the  gradients  from  internal  coordinates  to  the  Cartesian,  fixed  lab-frame 
coordinates.  The  result  is  a  set  of  6N  first-order  equations  that  must  be  solved  to 
propagate  the  trajectory. 
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There  are  numerous  integration  sehemes  that  have  been  devised  to  solve  sets  of 
differential  equations.  These  schemes  are  based  on  finite  steps  along  the  function  or 
system  of  functions  that  is  being  integrated.  These  finite  steps  are  used  to  extrapolate  to 
the  next  integration  point.  The  finite  nature  of  the  integration  process  means  that  a  stable 
integration  in  MD  is  determined  by  having  a  smooth  differentiable  potential,  since  the 
gradient  of  potential  is  used  in  Eq.  (10)  to  determine  how  the  momentum  changes. 
Changes  in  the  momenta  in  turn  determine  the  next  displacement  in  configuration  space 
at  which  the  potential  and  gradient  are  evaluated,  and  the  gradient  in  turn  affects  the 
momenta,  and  so  the  process  repeats.  Continuing  this  stepwise  process  is  known  as 
integrating  a  trajectory.  The  largest  cost  of  an  MD  simulation  is  typically  found  in  the 
evaluation  of  the  forces,  which  can  become  laborious  if  the  evaluation  is  done  by  finite 
difference  or  if  the  PES  is  complex,  as  in  direct  dynamics  or  interpolated  surfaces.  Even 
a  large  number  of  pairwise  interactions  to  be  evaluated  can  result  in  unfavorable  scaling 
of  the  computational  cost  of  a  problem.  Ultimately,  this  means  that  minimizing  the 
number  or  cost  of  gradient  evaluations  per  integration  step  leads  to  a  more  efficient 
simulation.  A  simple  means  of  decreasing  the  computational  cost  of  a  trajectory  is  to  use 
a  larger  integration  step  size,  so  that  fewer  integration  steps  are  required  to  simulate  some 
time  interval;  however,  but  care  must  be  taken  because  larger  step  sizes  can  increase  the 
effect  of  truncation  errors  . 

The  original  simulation  code  (GENDYN),  which  was  further  developed  in  this 
research,  had  the  4*  order  Runge-Kutta-Gill^^  (RKG4)  algorithm  as  the  standard 
integration  algorithm.  This  method  is  attractive  because  of  the  relatively  high  accuracy  of 
the  method,  0(h‘^),  and  is  self-starting  in  that  the  algorithm  does  not  require  multiple 
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forward  and  backward  integration  steps  to  calculate  the  next  step.  However,  this 
integrator  requires  four  gradient  evaluations  for  a  single  time  step,  so  direet  dynamies  or 
some  other  eostly  source  for  gradient  information  becomes  cost  prohibitive  with 
increased  integration  time  or  system  size. 

Numerieal  errors  in  integration  ean  easily  aeeumulate  and  lead  to  instability  in 
the  integration.  Pritehard"^'^  diseussed  the  effects  of  integration  step  size  and  finite 
preeision  on  the  reprodueibility  of  trajectories  on  various  maehines.  One  of  the  findings 
of  the  study  was  that  while  aeeeptable  agreement  of  ensemble  averages  aeross  the  various 
maehines  was  found,  results  of  individual  trajeetories  could  vary  substantially  due  to 
numerieal  errors.  Care  must  be  taken  in  considering  the  errors  of  finite  mathematics  that 
can  creep  into  a  simulation.  A  simple  illustration  of  the  issue  of  eumulative  error  ean  be 
seen  in  Fig.  3.  These  are  trajeetories  integrated  on  the  Henon-Heiles  potential  with  a  total 
energy  of  0.0125  using  two  different  integration  methods:  RKG4  and  veloeity  Verlet"^^ 
with  the  same  integration  step  size  of  0.1.  The  total  integration  time  was  1,000,000  time 
units.  All  of  the  trajeetories  were  started  at  the  same  point  in  phase  space.  On  a  shorter 
timescale,  frames  (a)  and  (b)  show  that  the  energy  drift  for  the  velocity  Verlet  and  RKG4 
integrators  is  very  small  for  the  entire  integration  time.  However,  as  the  RKG4  trajeetory 
is  integrated  for  a  mueh  longer  time,  there  is  a  progressive  deterioration  of  the  energy 
eonservation.  This  deterioration  is  also  refieeted  in  the  surfaces-of-seetion  plotted  in 
frames  (e-h)  for  the  long  time  trajeetories,  where  the  progressive  loss  of  energy 
eonservation  results  in  the  RKG4  trajectory  exploring  different  regions  of  phase  space. 
The  surfaees-of-seetion  were  collected  over  the  entire  integration  time.  Reduetion  of  the 
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integration  step  size  in  the  case  of  RKG4  by  an  order  of  magnitude  results  in  improved 
energy  conservation,  but  this  comes  with  the  additional  cost  of  more  force  evaluations. 
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Figure  4.  Plots  of  integration  accuracy  tests  for  the  velocity  Verlet  integrator  and  RKG 
with  identical  step  sizes  of  0.1:  (a)  Energy  drift  of  a  velocity  Verlet  trajectory,  (b)  same  as 
(a)  except  for  a  RKG  trajectory;  (c)  Energy  drift  of  a  velocity  Verlet  trajectory  plotted  in 
(a)  for  longer  timeframe,  (d)  same  as  (c)  except  for  a  RKG  trajectory;  (e)  Surfaces-of- 
section  for  the  X,  Px  subspace  from  a  velocity  Verlet  trajectory,  (f)  same  as  (e)  except  for 
RKG  trajectory;  (g)  Surfaces-of-section  for  the  Y,  Py  subspace  from  a  velocity  Verlet 
trajectory;  (h)  same  as  (g)  except  from  a  RKG  trajectory. 

The  potential  issue  of  loss  of  energy  conservation  does  not  mean  that  the  RKG 
integrator  is  not  useful,  since  the  order  of  accuracy  for  short  integrator  times  is  superior  to 
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that  of  the  Verlet  algorithm.  The  stability  of  the  Verlet  integration  is  due  to  it  being  part 
of  a  class  of  integrators  known  as  symplectic  integrators.  A  symplectic  integrator  is  one 
that  preserves  the  volume  of  phase  space  and  the  constants  of  motion  of  the  system."^^ 
Gray  et  tested  several  symplectic  integrators  and  concluded  that  there  is  not 

necessarily  one  integrator  that  is  the  best  for  all  applications. 

The  inefficient  scaling  of  the  RKG4  with  system  size  and  force  evaluations 
motivated  the  coding  of  two  additional  integrators  into  GENDYN;  one  integrator  being 
the  Adams-Moulton  Predictor-Corrector  algorithm"^^  and  the  other  the  velocity  Verlet. 
Each  of  these  integration  methods  was  coded  for  a  specific  use.  The  predictor-corrector 
was  coded  since  it  has  the  same  order  of  accuracy  as  RKG4,  which  is  important  for 
detailed  energy  transfer  studies,  but  only  requires  one  gradient  evaluation  per  integration 
step  as  opposed  to  four.  The  velocity  Verlet  was  coded  for  larger  systems  and  longer 
timescales  where  stable,  efficient  long  time  integration  becomes  important. 

E.  Sampling  and  Analysis 

The  majority  of  the  analysis  techniques  are  discussed  in  subsequent  chapters  in 
the  context  of  the  research  topic.  However,  there  are  a  couple  of  issues  that  are  worth 
discussing,  but  do  not  readily  fit  in  the  other  chapters  of  this  dissertation. 

Any  number  of  properties  can  be  calculated  from  a  classical  trajectory  simulation, 
but  there  is  an  important  caveat  that  must  considered  when  interpreting  simulation 
results.  Individual  trajectories  can  provide  mechanistic  details,  but  the  meaningful  result 
is  the  ensemble  average.  The  mean  of  the  ensemble  average  is  extremely  important  when 
considering  the  correspondence  between  classical  and  quantum  mechanics. Sewell 
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et  and  Guo  et  al.^^  both  studied  the  issue  of  zero-point  energy  flow  and  showed  that 
both  active  and  passive  constraints  on  zero  point  energy  result  in  aphysical  dynamics. 
Proton  tunneling"^^  can  also  be  an  important  quantum  effect  to  consider  when  studying  the 
dynamics  of  HO2  in  the  case  of  isomerization.  A  detailed  study  of  this  phenomenon  is 
beyond  the  scope  of  the  research  in  this  dissertation,  and  thus  such  a  study  is  outlined  as  a 
proposed  future  work  at  the  end  of  chapter  III. 

Error  estimation  in  trajectory  studies  can  be  done  with  two  different  approaches. 
The  first  approach  is  to  averaging  across  multiple  independent  groups  of  trajectories.  This 
approach  is  time  and  labor  intensive,  and  fortunately,  the  second  way  is  more  tractable. 
The  more  tractable  approach  is  known  as  bootstrap  estimation,  which  was  originally 
proposed  by  Efron.  The  bootstrap  procedure  consists  of  randomly  resampling  the  data 
and  generating  bootstrap  data  sets,  or  resampled  data  sets.  The  random  resampling  is 
done  through  use  of  a  pseudo-random  number  generator.  The  resampled  data  sets  are  then 
subject  to  the  same  analysis  as  the  parent  data  and  resampled  parameters  are  then  used  to 
calculate  statistics  that  can  be  used  to  describe  the  parent  data  set. 

We  applied  the  bootstrap  method  in  resampling  of  trajectory  indices  to  generate 
bootstrapped  ensembles  of  trajectories  for  analysis.  Our  application  of  the  bootstrap 
method  in  both  Chs.  Ill  and  IV  differed  from  the  original  bootstrap  method.  The  first 
difference  is  that  we  chose  to  have  the  bootstrapped  ensembles  of  trajectories  be 
composed  of  fewer  trajectories  than  the  parent  trajectory  ensemble.  The  traditional 
application  of  the  bootstrap  involves  having  the  number  of  points  in  each  resampled  data 
set  equal  to  the  number  of  points  in  the  parent  data  set.  It  is  well  known  that  the 
magnitude  of  the  statistical  error  of  an  average  calculated  in  a  trajectory  simulation  is  ~ 
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l/ ^^n  ,  where  n  is  the  number  of  trajectories  in  a  simulation.  This  means  that  bootstrapped 
data  composed  of  smaller  subsets  of  the  parent  data  should  provide  an  upper  estimate  to 
the  error  of  the  simulation. 

The  second  difference,  which  only  applies  to  Ch.  IV,  is  that  the  bootstrap  was 
modified  so  resampling  was  done  without  replacement,  to  prevent  duplicate  inclusion  of 
a  single-trajectory  index  in  an  individual,  bootstrapped  data  set.  Although,  it  is  possible 
that  a  trajectory  index  could  have  been  included  in  multiple,  independent  bootstrapped 
data  sets.  This  modification  to  the  resampling  was  done  because  a  trajectory  index  is  just 
a  placeholder  for  an  entire,  correlated  time  history.  However,  it  is  not  completely  clear  to 
us  if  resampling  of  smaller  data  subsets  without  replacement  is  equivalent  to  the  original 
bootstrap  method;  consequently,  we  treat  our  error  bars,  in  Ch  IV,  as  rough-upper 
estimate  of  the  variance  in  the  data. 
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III.  Dynamics  and  Unimolecular  Dissociation  of  Hydrogen  Peroxyl 
Radical 

A.  Introduction 
1,  Literature 

The  hydrogen  peroxyl  radical  HO2  is  important  in  free  radical  chain  propagation  and 
termination  in  combustion  reactions. Over  the  past  quarter  century  the  bimolecular 
reactions 


H  +  02^H0  +  0  (12) 

HO  +  O^H+02  (13) 

and  associated  unimolecular  dissociation  reactions 

HO2*  ^  H  +  O2,  (14) 

H02*^H0  +  0  (15) 


have  been  the  subjects  of  a  number  of  experimental  and  theoretical  studies.  An  important 
part  of  these  studies  has  been  potential  energy  surface  (PES)  development  using  ab  initio 
and  the  use  of  these  PESs  for  rate  calculations.  Brandao  et  al^'' 
reviewed  these  HO2  PESs  with  the  exception  of  the  most  recent  one  by  Ei  et  al.^^ 

We  have  used  classical  trajectories  to  study  the  fundamental  reaction  dynamics  of 
HO2  because  two  new  PESs  are  available,  both  based  on  accurate,  state-of-the-art 
electronic  structure  theory.  The  most  recent  PES,  by  Ei  et  al,  is  an  interpolating  moving 
least  squares  (IMLS)  polynomial  surface  that  constructed  from  ab  initio 
icMRCI+Q^^’^^’^°  results  with  CBS  extrapolation.^^  The  second  PES  is  the  third  member 
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of  a  progression  of  3-D  cubic-spline-fitted  surfaces  started  by  Xu,  Xie,  Zhang,  Lin,  and 
(XXZLG).  This  PES^^  is  based  on  =^1 8,000  icMRCI+Q  calculations,  but 
without  CBS  extrapolation.  It  best  represented  the  experimental  vibrational  spectroscopy 
of  the  PESs  discussed  in  the  2009  Brandao  et  review.  We  contrast  the  results 
calculated  using  these  recent  PESs  with  those  obtained  by  using  the  semi-empirical 
double  many-body  expansion  (DMBE-IV)  PES  developed  by  Pastrana  et  a/.^° 

The  DMBE-IV  PES  has  been  used  in  many  studies  of  the  reactions  dynamics  of  HO2. 
It  was  built  with  data  from  several  sources,  including  the  ab  initio  data  reported  by 
Melius  and  Blint  and  Walch  et  al.  ,  with  parameters  adjusted  to  match  experimentally 
determined  force  constants.  The  ab  initio  data  from  Walch  et  al.  for  the  O  +  OH 
asymptote  were  semi-empirically^"^  adjusted  for  electron  correlation.  The  method  used  to 
construct  the  DMBE-IV  PES  has  been  discussed  in  detail  elsewhere. Comparisons  of 
different  regions  of  the  PES  to  higher  quality  electronic  structure  calculations^ than 
used  in  its  construction  have  shown  significant  differences,  and  comparisons  to 
spectroscopic^^  and  kinetics^^’^^  experiments  show  some  of  the  limitations  of  this  surface 
in  accurately  describing  this  radical.  Several  classical  trajectory  studies^^’^*^  and  quantum 
dynamical  studies  ’  ’  of  microcanonical  unimolecular  decay  or  bimolecular  reaction 
dynamics  on  this  PES  have  shown  strong  mixing  among  the  modes,  and  reaction  kinetics 

77 

qualitatively  consistent  with  statistical  theories.  However,  a  trajectory  study  of  Miller 
using  the  Melius  and  Blint^^  PES  and  a  follow-up  trajectory  study  of  Miller  and  Garrett^"^ 
using  the  DMBE  IV  showed  that  the  HO2*  complexes  formed  from  H  +  O2  or  OH  +  O 
have  biased  decay  back  to  the  reactants.  This  result  implies  the  collision  complex  retains 
some  memory  of  its  origin.  Also,  Duchovic  and  Parker^^  inferred  from  a  classical 
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trajectory  study  of  H  +  O2  that  IVR  in  HO2*  is  slow  enough  to  cause  non-statistical 
unimolecular  decay,  an  observation  made  earlier  by  Lemon  and  Hase^^  based  on  an 
empirical  PES.  Marques  and  Varandas  examined  mode-specific  unimolecular 
dissociation  at  constant  energy  and  found  that  the  energy  decay  rate  constants  were  only 
weakly  dependent  on  the  mode  initially  excited.  However,  the  results  were  reported  at 
insufficient  resolution  to  compare  to  the  modest  effects  Miller  and  Garrett^"^  found.  The 
finding  of  these  studies  suggest  that  mode  mixing  is  strong  enough  on  the  DMBE-IV  PES 
to  create  an  irregular  spectrum  of  quantum  mechanical  metastable  states  and  a  single  time 
scale  of  unimolecular  decay,  but  not  strong  enough  to  completely  erase  all  memory  of  the 
initial  conditions  under  which  HO2*  is  formed. 

Studies  on  the  XXZEG  show  that  the  modes  of  HO2  are  more  quasi- 

periodic  and  decoupled  than  the  modes  predicted  using  the  DMBE-IV  PES.  The 
vibrational  spectrum  for  the  bound  states  on  the  XXZEG  calculated  by  Ein  et  alJ’’ 
showed  a  reduced  density  of  states  supported  on  the  surface,  with  the  distribution  of 
states  being  an  intermediate  case  of  regular  and  irregular  level  spacing.  They  found  that 
the  statistical  theory  rate  overestimated  the  average  quantum  mechanical  rate,  and  posited 
that  the  result  is  a  consequence  of  slow  IVR.  Xu  et  al.  found  regularities  in  the 
vibrational  spectrum  of  the  radical  even  at  energies  near  the  dissociation  threshold 
indicative  of  weak  modal  coupling,  which  they  suggested  would  strongly  affect 
unimolecular  dissociation  of  HO2.  They  assigned  pure  0-0  vibrational  overtones  up  to 
18  quanta  and  to  within  ~100  cm'^  of  the  dissociation  limit.  Pure  0-H  vibrational 
overtones  could  be  assigned  up  to  4  quanta,  but  at  higher  energies  isomerization  appears 
to  enhance  the  interaction  between  the  OH  stretching  and  OOH  bending  motions,  thus 
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ending  the  overtone  progression.  Xu  et  al.  also  suggested  that  resonanees  between  0-0 
and  0-H  could  also  play  a  role  in  energy  transfer^^*^^^^  between  these  two  modes.  The 
results  of  three  different  quantum  dynamics  studies  ’  on  XXZLG  for  the  bimolecular 
reactions  (la)  and  (lb)  show  significant  non-statistical  effects  in  the  differential  cross 
sections,  trace  the  contributions  from  individual  and  assignable  HO2*  resonances  to 
reaction  cross  sections  and  display  significant  differences  with  quantum  scattering  results 
on  DMBE-IV.  Quasiclassical  trajectory  calculations  on  the  XXZLG  PES  have  been 
carried  out  for  reactions  (la)  and  (lb)  by  Lendvay  et  al.  and  show  that  the  lifetimes  of 
the  collision  complex  formed  in  reaction  (la)  display  a  bi-exponential  distribution.  Eor 
reaction  (lb),  Jorfi  et  al.  found  evidence  of  slow  energy  transfer  between  the  0-H  and 
0-0  bonds,  resulting  in  a  non-statistically  large  probability  of  HO2*  dissociating  back  to 
the  H  +  O2  reactants.  All  of  these  studies  suggest  that  XXZLG  dynamics  has  much 
weaker  mode  coupling  than  DMBE-IV  dynamics,  and  the  result  is  much  more  salient 
non-statistical  behavior  in  bimolecular  reactions.  Despite  the  high-quality  ab  initio 
foundation  of  XXZLG,  exact  quantum  dynamical  rate  constants  on  XXZLG  for  the  O  + 
OH  reaction  (lb)  agree  with  experiment  within  a  factor  of  two.  Likely  sources  of  the 
discrepancy  include  experimental  error,  limitations  of  the  PES,  and  the  possible 
relevance  of  multistate  non-adiabatic  dynamics. 

The  possible  involvement  of  electronically  excited  states  and  the  need  for  more 
accurate  electronic  structure  methods  led  to  the  calculations  of  new  sets  of  PESs  for  the 
HO2  ground-  and  first-excited  states  by  Li  et  al.^^  Among  those  PESs  is  the  ground  state 
IMLS  PES  that  is  further  studied  in  this  paper.  The  electronic  structure  theory  used  to 
generate  the  IMLS  PES  includes  extrapolation  to  estimate  the  CBS  limit.  While  also 
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interpolative,  IMLS  fitting  is  a  different  kind  of  representation  than  the  splines  used  in 
XXZLG  and  is  eapable  of  high  accuracy  with  relatively  few  points.  However,  while 
containing  the  reaction  path  for  dissociation  to  H  +  O2,  the  current  IMLS  PES  is  restricted 
to  coordinate  ranges  that  do  not  include  the  higher  energy  O  +  OH  asymptote. 

This  brief  survey  of  the  literature  shows  that  earlier  work  placed  some  emphasis  on 
the  extent  of  mode  coupling  in  HO2*.  Despite  that  emphasis,  lower-energy  processes, 
such  as  IVR  and  isomerization,  which  can  be  sensitive  to  mode  coupling,  have  not  been 
extensively  studied  on  any  of  the  three  surfaces.  If  lower-energy  processes  display 
different  results  on  the  different  PESs,  perhaps  it  would  be  possible  to  identify  specific 
PES  features  that  underpin  the  dynamical  differences.  Eor  the  higher-energy  process  of 
unimolecular  decay,  mode-specificity  has  only  been  studied  on  the  DMBE-IV  PES  where 
mode  coupling  is  thought  to  be  quite  strong  and  only  modest  mode-specific  effects  have 
been  identified.  Based  on  previous  IVR  studies,^"^*^^^’'*’^  mode-specific  effects  in 
unimolecular  decay  are  expected  to  be  more  significant  for  the  XXZEG  PES  (due 
apparently  to  weaker  mode  coupling).  No  dynamics  studies  have  been  carried  out  using 
the  IMES  PES. 

In  this  study  we  investigate,  using  classical  trajectories,  the  intramolecular  dynamics 
and  isomerization  on  all  three  PESs  and  the  unimolecular  dissociation  of  HO2  ^  H  +  O2 
on  the  XXZEG  and  DMBE-IV.  This  chapter  is  organized  as  follows:  Descriptions  of  the 
characteristics  of  the  PESs  are  given  in  Sec.  A.2,  the  computational  methods  are 
described  in  Sec.  B,  results  and  discussion  are  given  in  Sec.  C,  and  the  conclusions  are 
given  in  Sec.  D.  Potential  future  studies  are  proposed  in  Sec.  E. 
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2.  Potential  Energy  Surfaces 

a.  The  XXZLG  PES 

Four  different  spline -based  versions  of  the  XXZLG  PES  have  been  produeed  by  Xie 
and  eoworkers.  The  first  two^^’^"^  are  based  on  -15,000  ab  initio  data  points  computed  at 
the  level  with  the  valence  electrons  correlated  in  the  icMRCI+Q/AVQZ  calculations.  The 
active  space  for  the  CASSCF  reference  (9e,  7o)  constrains  the  oxygen  I5  and  2s  orbitals 
to  be  doubly  occupied  and  includes  two  A"  states  with  equal  weights.  Later  in  2007  the 
data  set  was  augmented  by  =3,000  points  for  large  0-0  separations  (as  encountered  in 
the  O  +  OH  channel),  leading  to  the  third  version  of  the  XXZLG  PES.^^  Finally,  in  2010 
as  part  of  a  study  of  the  lowest  A'  excited  state,  a  further  update  to  the  ground  state  data 
was  also  reported.  The  new  ab  initio  data  were  computed  using  a  full  valence  active 
space  (13e,  9o)  in  the  two-state  CASSCF  reference  and  the  core  electrons  were  also 
correlated  in  the  subsequent  icMRCI+Q  calculations.  Including  core  electron  correlation 
and  the  full-valence  active  space  was  reported  to  improve  the  thermochemistry  slightly 
with  respect  to  experiment.  For  the  ground  state,  these  improved  ab  initio  calculations 
were  carried  out  at  =4,000  points  to  obtain  their  difference  with  respect  to  the  previous 
XXZLG  PES.  This  difference  was  fit  with  a  spline  which  when  added  to  the  previous 
XXZLG  PES  constitutes  the  fourth  version  of  the  XXZLG  PES.  The  version  of  the 
XXZLG  PES  used  in  this  study  is  the  third  version^^  developed  in  2007  from  =18,000  ab 
initio  points. 

An  examination  of  the  isomerization  saddle  point  identified  a  minor  seam  in  the 
XXZLG  PES  that  results  in  a  discontinuous  force  for  trajectories  traversing  the  0-0 
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perpendicular  symmetry  plane.  This  seam  occurs  because  XXZLG  uses  only  the  shorter 
of  the  two  0-H  distances.  The  use  of  the  shorter  0-H  distance  guarantees  a  continuous 
value  across  the  symmetry  plane  but  not  a  continuous  gradient.  We  smoothed  the  seam 
using  symmetrization  ’  of  the  PES  in  the  region  near  the  symmetry  plane.  This 
smoothed  surface  has  the  form 

-^(H-0,0J5(A)  +  E(H-0,0J(1-5(A))  (16) 

where  ^(A)  is  a  smoothly  varying  second-order  switch  of  the  following  form 

S(A)  =  min[l,  max[0,  1-  (lOg(A)'  -  15g(A)"  +  6g(A)')]]  (17) 

where 

g(A)  =  (A-Ai)(A-A,)-i  (18) 

where  the  progress  variable  A=7?H-Oa  -  f?H-ob  with  limits  A^  =  0.2  A  and  A/  =  -0.2  A.  The 
value,  gradient,  and  Hessian  of  ^(A)  are  zero  at  both  limits,  making  the  V^ym  patch 
smoothly  join  the  unmodified  XXZLG  PES  outside  of  -0.2  A  from  the  0-0 
perpendicular  bisector  plane.  Note  also  that  precisely  at  the  symmetry  plane,  the  patched 
surface  is  exactly  equal  to  the  original  XXZLG  PES.  Reproducing  the  value  of  the 
original  XXZLG  surface  at  the  symmetry  plane  means  that  the  isomerization  saddle  point 
energy,  geometry,  and  vibrational  frequencies  are  exactly  the  same  on  the  original  and 
corrected  XXZLG  PES  except  for  the  imaginary  frequency,  which  is  ill-defined  on  the 
original  XXZLG  PES.  Additional  technical  details  of  the  XXZLG  PES  are  discussed  in 
the  appendices. 
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b.  The  IMLS  PES 

Li  et  al.^^  reported  a  ground-state  IMLS  PES  for  HO2  based  on  icMRCI+Q  data  using 
CBS  extrapolation  and  18-referenoe-states  with  dynamic  weighting.  We  used  the  Cs 
point  group  and  built  the  set  of  molecular  states  from  separate  ground  state  atoms, 
H02('A")  ^  O('Pg)  +  O('Pg)  +  R{%).  (19) 

Resolving  the  atomic  states  into  the  molecular  Cs  point  group  using  tables  from 
Herzberg  and  combining  into  molecular  states  yields  a  total  of  18  doublet  states  (10  A', 
8  A")  which  are  asymptotically  degenerate  for  separated  atoms.  There  are  also  sets  of 
higher  spin  quartet  and  sextet  states  not  considered  here.  The  18  doublet  states  were 
included  in  dynamically  weighted  (DW)  state-averaged  (SA)  CASSCF  calculations^^’^^ 
using  a  full-valence  active  space.  Maximum  weight  was  focused  on  the  ground  A"  state, 
and  weights  for  the  other  states  were  determined  by  the  DW  scheme  with  B  =  2.0  eV.^*^ 
The  DW  scheme  correctly  reflects  the  differing  degeneracies  which  arise  in  different 
regions  of  the  PES,  including  Renner-Teller  degeneracy  with  the  lowest  A'  state 
observed  at  linear  geometries.  The  DW  multistate  scheme  also  promotes  robust 
convergence  that  is  critically  important  for  automated  interpolative  methods  (such  as 
IMES)  in  which  new  data  is  added  automatically.  Using  the  DW  18-state  description,  no 
convergence  problems  such  as  those  noted  in  the  2-state  fixed  weight  calculations 

c-j 

reported  by  Xie  et  al.  were  encountered.  The  18-state  DW-SA-CASSCF  calculations 
were  used  as  the  reference  for  subsequent  icMRCI+Q  calculations. All  electrons  were 

-5 

correlated  in  the  icMRCI  calculations,  and  the  CBS  limit  was  estimated  using  the  /' 
formula^^  and  Dunning’s  aug-cc-pVTZ  and  aug-cc-pVQZ  bases. 
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An  automated  PES  generation  seheme^^*^'’^”^^’^^  was  used  to  construct  the  IMLS  PES. 
A  total  of  1 ,609  symmetry-unique  points  was  generated  by  the  algorithm  starting  with  a 
seed  grid  of  200  points.  Jacobi  coordinates  were  used  with  ranges  of  Rqo  =  [1-0,  2.1]  A 
and  Rh-oo  =  [0.5,  2.65]  A.  Data  were  computed  for  y  <  nil  with  the  symmetry  partner 
(7i-y)  added  at  no  cost.  Explicit  inclusion  of  the  symmetry  partner  ensures  correct 
exchange  symmetry  and  provides  smooth  behavior  across  the  transition  between  the  two 
wells.  Energy  was  restricted  to  a  maximum  of  70  kcal/mol  above  the  global  minimum. 
The  algorithm  for  generating  the  PES  as  well  as  the  form  of  the  interpolative  weight 
function  and  related  numerical  parameters  are  as  described  in  the  supporting  information 
of  Ref.  28(b).  The  automatic  PES  generation  was  terminated  with  1,609  points  when  the 
estimated  fitting  error  reached  1  cm'\  The  computed  ab  initio  value  of  the  isomerization 
saddle  barrier  height  differs  by  2  cm'^  from  the  IMES  PES  value  while  the  geometries 
differ  by  less  than  0.0001  A,  consistent  with  the  estimated  quality  of  the  fit.  It  has  been 
shown  in  previous  studies  that  once  an  IMES  PES  is  fit  to  negligible  fitting  error 
(~lcm'*),  very  high  fidelity  to  the  structural  parameters  and  vibrational  levels  associated 
with  that  ab  initio  method  is  achieved.^^ 

c.  Comparisons  of  the  IMLS,  XXZLG,  and  DMBE-IV  PESs 

Eor  all  three  PESs,  the  dissociation  energy  (both  Do  and  De)  and  the  isomerization 
saddle  point  energy  are  listed  in  Table  I  along  with  the  measured  dissociation  energy. 
The  best  experimental  estimate  for  Do(H-OO)  is  48.01±0.04  kcal/mol. The  value  of 
Do(H-OO)  is  0.21  kcal/mol  lower  for  the  IMLS  PES,^^  0.82  kcal/mol  lower  for  XXZLG, 
and  0.73  kcal/mol  higher  for  DMBE-IV.  Typically  ab  initio  methods  cannot  capture  all 
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the  correlation  energy  of  a  bond  and  thus  underestimate  dissociation  energies.  The 
remaining  -0.2  kcal/mol  discrepancy  in  the  IMLS  result  is  assigned  to  error  in  the  CBS 
procedure  as  well  as  the  effect  of  higher-order  correlation  not  captured  by  the  icMRCI 
method.  The  classical  isomerization  barrier  is  highest  on  the  DMBE-IV  PES,  with  IMES 
0.798  kcal/mol  lower  and  XXZEG  2.01  kcal/mol  lower  than  IMES.  The  XXZEG  value  in 
Table  I  is  the  value  for  the  fitted  PES.  As  discussed  above,  when  the  electronic  structure 
method  that  XXZEG  is  based  on  is  used  to  directly  calculate  the  barrier,  the  result  is  1 .3 
kcakmol  higher  than  the  value  in  the  table.  This  barrier  is  -0.7  kcal/mol  lower  than  the 
IMES  barrier  and  consistent  with  the  differences  found  for  dissociation. 

The  entries  in  the  lower  part  of  Table  I  are  the  geometries  of  the  equilibrium  and 
isomerization  saddle  point  structures.  The  experimental  HO2  equilibrium  geometry  is 
essentially  exactly  given  by  the  DMBE  PES  (which  was  adjusted  to  do  this)  and  closely 
approximated  by  the  XXZEG  and  IMES  PESs.  While  the  XXZEG  and  IMES  PESs  have 
very  similar  isomerization  saddle  point  geometries,  the  DMBE  PES  has  a  noticeably 
more  elongated  structure;  that  is,  longer  Rqq  and  7?oh  distances. 
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Table  I  Dissociation  and  isomerization  energies  and  equilibrium  and  isomerization 
saddle  point  geometries  for  the  DMBE-IV,  XXZLG,  and  IMLS  PESs;  and  available 
experimental  results 


DMBE-IV 

XXZLG 

IMLS 

Experiment 

Energetics  (kcal/mol) 

47.19’’’" 

47.80’’’"’‘’ 

Dissociation  to  H+O2  Do 

48.74^^ 

48.01±0.04" 

De 

54.85 

53.68 

54.36 

Isomerization  Barrier 

40.72 

36.92 

38.93 

Geometries 

Equilibrium 

0.9705’ 

r(OH)  A 

0.9708 

0.9708 

0.9689 

r(00)  A 

1.3305 

1.3341 

1.327 

1.330’ 

(OOH)" 

104.29 

104.12 

104.39 

104.29’ 

Isomerization  Saddle  Point 

r(OH)  A 

1.202 

1.174 

1.162 

r(00)  A 

1.484 

1.424 

1.420 

(OOH)" 

51.88 

52.66 

52.34 

^  Mandelshtam  et  al.,  Ref  66. 

’’  The  O2  harmonic  ZPE  required  for  DO  is  787  om'\  There  is  an  inconsequential 
anharmonic  contribution  to  the  O2  zero  point  energy;  experimentally  determined  WeXg 
is  ~12  em'*  [P.  H.  Krupenie,  J.  Chem.  Ref.  Data  1,  423  (1972)]. 

The  HO2  ZPE  used  in  ealculating  for  Do  has  an  anharmonic  value  of  3054.6  cm'^  (for 
XXZEG)  and  3083.2  cm'^  (for  IMES)  calculated  as  described  in  the  text. 

The  fitted  range  of  the  IMES  PES  only  extends  to  Rh-oo  =  2.65  A,  so  the  asymptotic 
energy  was  evaluated  from  an  eleetronie  structure  caleulation  with  Rh-oo  =  25.0  A. 
The  difference  in  energy  between  25.0  A  and  2.65  A  is  0.1 1  kcal/mol. 

®  Ruscic,  Ref  97. 

^  Eubie,  Ref  63(a). 

Table  II  lists  the  ealculated  lower  vibrational  levels  on  the  three  PESs  and  the 
experimental  values.  Eor  the  IMES  and  XXZEG  PESs,  the  entries  in  Table  II  and  the 
HO2  zero  point  energy  (ZPE),  referred  to  in  Table  I,  are  from  a  potential-optimized 
discrete  variable  representation  (PODVR)  variational  ealculation.^^  The  details  of  those 
calculations  are  as  reported  in  Ref  56  for  the  IMLS  PES.  Vibrational  levels  have  been 
reported^^’^"^  for  the  two  earliest  versions  of  XXZLG  PESs  but  not  for  the  later  version 
used  here.  The  XXZLG  and  IMLS  PESs  predict  vibrational  spectra  in  excellent 
agreement  with  experiment.  Eor  the  seven  levels  listed  in  Table  II,  the  root-mean-square 
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(rms)  error  with  respect  to  experiment  is  ~9  cm'^  for  the  XXZLG  PES  and  ~10  cm"'  for 
the  IMLS  PES.  In  contrast,  the  rms  error  for  the  DMBE-IV  PES  is  an  order  of  magnitude 
larger. 

We  have  found  previously  that  a  litany  of  small  corrections  such  as  core  correlation, 
CBS  extrapolation,  spin-orbit  or  relativistic  corrections  can  offset  each  other,  and  unless 
everything  is  included,  improving  one  aspect  of  the  electronic  structure  calculation  might 
not  improve  agreement  with  experiment  (if  the  lower-level  calculations  benefit  from 
cancelation  of  errors).  Eurther  improvements  of  the  XXZEG  and  IMES  PESs  with  high- 
level  electronic  structure  theory  results  will  require  addressing  the  effects  of  these  many 
small  corrections. 


Table  II.  Comparison  to  experiment  of  the  calculated  vibrational  levels  in  cm'^  for  the 
DMBE-IV,  XXZEG,  and  IMLS  PESs. 


DMBE-IV" 

XXZLG 

IMLS 

Experiment 

H02 

(001) 

1065.5  (  32.1)” 

1090.1  (  7.5) 

1104.7  (  -7.1) 

1097.63‘^ 

(010) 

1296.4  (  95.4) 

1388.4  (  3.4) 

1395.9  (  -4.1) 

1391.75" 

(100) 

3333.7  (102.5) 

3442.6  (  -6.4) 

3448.9  (-12.7) 

3436.20‘* 

(200) 

6492.4  (158.8) 

6671.0  (-19.8) 

6672.8  (-21.6) 

6651.19" 

DO2 

(001) 

1015.0  (  5.2) 

1023.4  (  -3.2) 

1020.16” 

(010) 

1116.5  (  5.0) 

1122.3  (  -0.8) 

1121.47” 

(100) 

2552.7  (  -3.5) 

2551.7  (  -2.5) 

2549.22® 

Brandao  et  al.,  Ref  57. 

’’  The  numbers  in  parentheses  are  the  differences  in  the  experimental  and  calculated 


values. 

°  Burkerholder,  63(e) 
Yamada,  63(d) 

®  DeSain,  63(f) 

^  Uehara,  63  (b) 

®  Lubic,  63(a) 


The  contour  plots  of  the  three  PESs  in  Eig.  5  highlight  the  similarity  in  global  shapes 


of  the  IMLS  and  XXZLG  PESs  and  how  they  differ  from  the  DMBE-IV  PES.  In  these 
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plots,  the  Jacobi  coordinates  i?oo-H  and  4  0oo-h  vary  while  the  0-0  bond  distance  is  fixed 

at  the  value  in  the  HO2  equilibrium  geometry  of  each  respective  PES.  The  plots  of  the 
XXZLG  and  DMBE-IV  PESs  are  matched  side  by  side  along  the  0-0  perpendicular 
bisector  symmetry  line  with  the  plot  of  the  IMLS  surface  in  order  to  provide  direct 
comparison  of  the  surfaces.  The  two  symmetrically  related  equilibrium  positions  can  be 
seen  with  the  isomerization  barrier  separating  them.  At  the  top  of  each  plot,  the  two 
symmetrically  related  dissociation  channels  from  each  isomer  are  readily  seen.  At  the 
bottom  of  each  plot,  there  are  two  symmetrically  related  barriers  to  linearity  in  the  H-0-0 
bending  motion. 


Figure  5  Contour  plots  of  the  IMES  (a),  XXZEG  (b),  and  DMBE-IV  (c)  surfaces, 
with  respect  to  ROH  and  the  center  of  mass  angle  of  H  with  the  0-0  bond.  The  0-0 
bond  distance  is  fixed  at  the  HO2  equilibrium  geometry  of  each  PES.  Contour  spacing  is 
10  kcal/mol. 
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On  the  scale  of  Fig.  5,  the  IMLS  and  XXZLG  PESs  in  the  top  panel  appear  to  be  very 
similar;  however,  there  is  a  slight  mismatch  (~1  kcal/mol)  at  the  isomerization  symmetry 
line  and  at  the  barrier  to  linearity.  The  width  of  the  dissociation  channel  illustrated  by  the 
closeness  of  the  red  features  in  the  plot  is  probably  slightly  narrower  in  the  IMLS  PES. 
These  results  are  consistent  with  the  energy  differences  in  Table  I.  The  bottom  panel  of 
Eig.  5  shows  that  the  IMLS  and  DMBE-IV  PESs  are  clearly  different,  with  the  DMBE-IV 
having  a  more  curved  reaction  path  between  the  isomerization  barrier  and  the  equilibrium 
point,  a  considerably  lower  barrier  to  linearity,  and  a  narrower  dissociation  channel.  Yang 
and  Klippenstein  commented  on  the  lack  of  sufficient  ab  initio  input  into  DMBE-IV  for 
the  bending  component  of  the  dissociation  channel  and  the  possible  implications  for 
reliable  rate  constant  calculations. 

The  manifestations  in  the  dynamics  of  the  similarities  and  differences  in  the  three 
PESs  are  clearly  illustrated  with  SOS  plots. The  SOS  plots  shown  in  Eig.  6  were 
calculated  on  each  of  the  three  PESs  with  groups  of  50  trajectories  with  energies  of  10.0, 
20.0,  30.0,  and  40.0  kcal/mol  (trajectory  details  are  given  in  Sec.  B.I.).  The  2N-2 
mapping,'®®  which  defines  the  4D  subspace,  was  done  by  freezing  the  0-0  distance  at  the 
respective  equilibrium  value  for  each  PES  (as  in  Eig.  5).  Xu,  et  al.  found  evidence  of 
regularity  in  the  motion  of  the  0-0  mode  up  to  dissociation,  but  also  found  signs  of 
irregularity  in  the  motion  of  the  modes  associated  with  the  hydrogen  atom.  This  irregular 
motion,  which  could  be  associated  with  mode  mixing,  is  what  we  would  like  to 
understand.  Surfaces-of-section,  shown  in  Eig.  6,  makes  the  visualization  of  the 
accessible  phase  space  and  identification  of  regularities  in  trajectories  straightforward  and 
allow  for  direct  comparison  for  reduced  dimensionalities  of  the  dynamics.  We  defined 
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two  planes  in  Cartesian  spaee  using  the  equilibrium  coordinates  of  the  hydrogen  atom 
{xmin,ymin)-  Flgurc  6  shows  the  record  of  the  coordinate  x  and  conjugate  momentum 
each  time  a  trajectory  of  the  hydrogen  atom  passes  through  the  plane  defined  by  ymin-  An 


analogous  figure  (not  included  here)  for  y  and  Py  for  passage  through  the  plane  defined 


by  Xmin  shows  similar  results. 


Figure  6  Poincare  surfaces-of-section  for  x  and  Px  subspace  (see  text  for  details).  From 
the  top,  rows  are  results  for  energies  10,  20,  30,  and  40  kcal/mol.  From  the  left,  columns 
are  for  the  IMLS,  XXZLG,  and  DMBE-IV  PESs 

The  results  in  Eig.  6  for  the  IMES  and  XXZEG  PESs  are  similar  over  the  whole 
energy  range.  Except  at  the  highest  energy  of  40  kcal/mol,  both  show  a  regular 
phase-space  structure.  At  40  kcahmol  there  is  increasing  irregularity  in  the  phase  space. 
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but  the  major  regular  features  of  the  SOS  from  the  two  PESs  still  mateh.  In  eontrast,  the 
phase  space  of  the  DMBE-IV  is  largely  irregular  except  at  the  lowest  energy  of  10 
kcahmol,  which  is  just  above  the  ZPE.  The  features  of  the  regular  phase-space  structure 
do  not  match  well  with  those  of  IMLS  and  XXZEG,  consistent  with  noticeable 
differences  in  the  vibrational  frequencies  as  seen  in  Table  IT  The  highest  total  energy  in 
Pig.  6  remains  below  the  isomerization  barrier  on  all  three  PESs  for  the  fixed  0-0 
distance  of  the  SOS  calculations.  (Note  that  the  isomerization  barrier  heights  in  Table  I 
are  for  the  0-0  distance  relaxed.)  Related  results  to  those  shown  in  Pig.  6  have  been 
reported  by  Seidel  et  al}^^  for  the  DMBE-IV  PES  and  by  Barnes  and  Kellman''^^  for 
earlier  versions  of  the  XXZEG  PES  than  the  one  used  here. 

These  results  on  energetics,  geometry,  frequencies,  shape,  and  phase-space  structure 
emphasize  the  similarities  between  the  IMPS  and  XXZEG  PESs  and  their  common 
differences  with  the  pioneering  DMBE-IV  PES.  The  most  easily  quantified  difference  in 
the  IMLS  and  XXZEG  PESs  is  in  the  isomerization  barrier  heights,  which  is  probably 
mostly  due  to  a  local  sparsity  of  ab  initio  data  used  to  fit  the  XXZEG  PES  (see 
Appendices)  and  only  partially  due  to  the  different  accuracies  of  the  underlying  electronic 
structure  theory  methods.  At  greater  levels  of  detail,  there  are  other  differences  between 
the  IMLS  and  XXZEG  PESs  that  have  implications  for  dynamics  and  are  explored  in 
Sec.  C. 
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B.  Methods 

1.  Trajectory  Calculations 

All  of  the  simulations  were  performed  with  the  classical  trajectory  code  GenDyn. 
Groups  of  50  trajectories  were  calculated  for  the  SOS  plots  and  groups  of  1,000 
trajectories  were  calculated  to  determine  all  rate  constants.  The  trajectories  were 
integrated  using  a  fixed  step  size  5*V6*’’-order  Adams-Moulton  predictor-corrector 
integrator"^^  initiated  with  a  fourth-order  Runge-Kutta-Gill  routine.  The  SOS  trajectories 
were  integrated  for  5  ps.  The  integration  time  for  the  IVR  studies  was  10  ps.  The 
trajectories  for  studying  reactive  processes  were  integrated  until  either  a  reaction  criterion 
was  satisfied,  or  the  maximum  integration  time  of  10  ps  was  reached.  Hessian  matrices 

103 

for  normal  mode  and  IVR  analyses  were  calculated  by  fourth-order  finite  difference. 
The  angular  momentum  was  set  to  zero  in  all  simulations. 

Statistical  initial  conditions  were  prepared  for  total  energies  of  44,  46,  48,  50,  and  52 
kcal/mol  for  isomerization  and  ~60.0  kcal/mol  for  dissociation  using  the  ejjicient 
microcanonical  sampling  algorithm  as  implemented  in  GENDYN.  All  three 

atoms  were  moved  during  a  Markov  walk  started  from  the  equilibrium  geometry  with  a 
250,000-step  warm-up  walk,  followed  by  100,000-step  walks  between  trajectories.  The 
maximum  atomic  displacement  of  each  atom  was  0.2  A,  and  the  acceptance/rejection 
ratio  of  each  walk  was  in  the  range  of  0.4  to  0.6.  Statistical  initial  conditions  were 
selected  for  total  energies  of  10,  20,  30  and  40  kcal/mol  for  SOS  plots.  Only  the 
hydrogen  atom  moved  in  the  Markov  walk  with  a  25,000-step  warm-up  with  subsequent 
walks  of  10,000  steps.  The  Markov  walk  assigned  only  the  x  andy  components  of  the 
momentum  of  the  hydrogen  atom.  The  momenta  of  the  oxygen  atoms  and  the 
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momentum  in  the  z  direetion  of  the  hydrogen  atom  were  fixed  at  zero,  restrieting  the 
phase  spaee  to  four  dimensions. 

Quasielassieal  initial  eonditions^"^  were  used  for  initial  exeitation  of  one  normal  mode 
while  the  other  normal  modes  were  assigned  ZPE.  At  high  energies  signifieant  regions  of 
phase  space  may  not  be  well  described  by  the  normal  modes,  which  can  cause  small 
transient,  minor  effects  in  the  initial  trajectory  behavior.  However  we  find  that  these 
minor  effects  do  not  influence  the  conclusions  of  this  study.  The  mode-specific  initial 
excitations  are  listed  in  Table  III,  where  the  normal  modes  are  identified  by  their  major 
component  in  valence  coordinates. 

Table  III.  Initial  normal  mode  excitation  energies  (and  corresponding  fractional 
quantum  numbers)  and  ZPE  energies.  Excitation  quantum  numbers  are  in  parentheses. 


DMBE-IV 

IMES 

XXZEG 

Total  Energy  (kcal/mol) 

E  =  ~35.5  (below  isomerization) 

0-0 

28.6  (8.6) 

28.1 

(8.2) 

28.4  (8.4) 

0-0-H 

28.9  (7.0) 

28.9 

(6.5) 

29.0  (6.6) 

0-H 

31.9  (2.7) 

31.6 

(2.5) 

31.5  (2.5) 

E  =  =47.3  (above  isomerization) 

0-0 

40.9(12.5) 

40.1  (11.8) 

39.9(12.0) 

0-0-H 

40.6(10.0) 

40.5 

(9.3) 

40.5  (9.4) 

0-H 

43.8  (3.9) 

43.7 

(3.7) 

43.6  (3.7) 

E  =  =59.9  (above  dissociation) 

0-0 

51.88  (16) 

50.34(15.26) 

0-0-H 

52.19  (13) 

50.72(11.91) 

0-H 

54.75  (5) 

53.45  (4.593) 

Harmonic  ZPE 

0-0 

1.57 

1.63 

1.60 

0-0-H 

1.93 

2.07 

2.04 

0-H 

4.98 

5.27 

5.25 

The  normal  modes  are  identified  by  the  major  internal  coordinate  component. 


Eor  IVR  and  isomerization  studies,  fractional  quantum  numbers  were  used  to  give  the 
same  total  energy  to  within  ~1%  for  each  of  the  PESs.  We  have  repeated  the  Marque  and 
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70 

Varandas  trajectory  study  of  unimolecular  dissociation  on  the  DMBE-IV  PES.  Eor  a 
comparable  calculation  using  the  XXZEG  PES,  we  assigned  ZPE  to  the  unexcited  modes 
and  used  fractional  quantum  numbers  (given  in  parentheses  in  Table  III)  for  the  excited 
modes  to  obtain  the  same  total  energy  above  De  for  the  two  PESs.  The  difference  in  the 
average  total  energies  for  the  DMBE-IV  and  XXZEG  PESs  corresponds  to  the  difference 
in  De  (see  Table  I).  The  approximate  total  energy  listed  in  Table  III  is  the  mean  of  these 
two  average  energies. 

2.  Analysis  of  Trajectories 

Studying  the  energy  relaxation  of  an  excited  vibrational  mode  requires  a  means  of 
measuring  energy  localization,  which  can  only  be  done  approximately  due  to  ro- 
vibrational  coupling.  The  method  selected  for  monitoring  IVR  involves  decomposing  the 
kinetic  energy  by  standard  transformations  into  the  three  normal  mode  kinetic  energies 
and  then  assuming  that  the  fraction  of  total  kinetic  energy  in  a  mode  is  equivalent  to  the 
fraction  of  total  energy  in  a  mode.  In  this  approach,  the  inherent  ambiguity  in  assigning 
potential  energy  to  different  modes  is  avoided.  Since  the  normal  mode  eigenvectors 
mathematically  span  the  space  of  the  three  internal  coordinates,  the  kinetic  energy 
decomposition  is,  in  principle,  complete.  By  this  means,  the  energy  E/t)  in  the  mode 
is  calculated  from  the  trajectories.  To  us  help  interpret  Ei(t),  we  developed  an  analytical 
kinetics  model  of  energy  transfer  in  a  three  mode  system  that  is  based  on  two  principles: 
(1)  the  change  in  the  energy  of  a  mode  is  directly  proportional  to  the  amount  of  energy 
the  mode  currently  has  and  (2)  the  energy  of  each  mode  asymptotically  relaxes  to  the 
same  value,  i.e.,  the  total  energy  divided  by  the  number  of  modes.  As  discussed  in  the 


42 


59 


Appendix  II,  these  two  principles  lead  to  an  explicit  expression  for  Ei(t)  in  terms  of  two 
effective  IVR  rate  constants  kj  and  one  constant  c  that  affects  the  amplitude  of  the  two 
decay  processes.  These  expressions  will  be  compared  to  the  classical  trajectory  results. 

The  two  bond-breaking  processes  of  isomerization  and  unimolecular  dissociation  are 
defined  as  follows.  The  time  of  isomerization  was  defined  to  be  when  the  hydrogen 
atom  passed  through  the  symmetry  plane,  bisecting  the  0-0  bond,  provided  that  the 
shorter  0-H  stretching  motion  underwent  two  turning  points.  This  arbitrary  definition 
excludes  rapidly  recrossing  trajectories.  For  dissociation,  the  lifetime  was  taken  to  be  the 
time  of  the  last  inner  turning  point  for  the  shorter  OH  bond  distance  before  the  distance 
between  the  center  of  mass  of  O2  and  H  exceeded  5  A.  This  inner  turning  point  definition 
eliminates  the  effect  of  the  relative  velocity  of  escaping  H-atom  on  the  dissociation 
lifetime.  By  the  turning  point  definition  of  a  lifetime,  any  trajectory  that  dissociates  with 
no  inner  turning  point  must  be  back  integrated  in  time  to  find  its  last  inner  turning  point 
in  negative  time.  We  refer  to  these  rapidly  dissociating  trajectories  as  “ballistic” 
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trajectories. 

The  first-order  rate  constants  for  isomerization  and  dissociation  were  calculated  by 
least-squares  fitting  of  the  unisomerized  or  undissociated  trajectory  populations 
(determined  from  the  lifetimes)  to  the  bi-exponential  form: 

E{t)  =  A  exp  (“AO  +  B  exp  {—k^t)  (20) 

where  A  and  B  are  prefactors  and  ki  and  k2  are  rate  constants.  The  prefactors  were 
determined  by  linear  least  squares  fitting  and  the  rate  constants  were  determined  by  a  grid 
search.  In  many  cases,  the  second  term  in  Eq.  (20)  was  found  to  be  negligible.  For  some 
mode-specific  excitations,  for  an  initial  period  of  time  there  is  either  no  decay  rate  or  an 
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accelerating  deeay  rate.  The  reasons  for  such  initial  regions  are  varied^"^  and  will  be 
discussed  as  they  appear  in  the  next  section.  We  will  refer  to  this  as  a  latency  effect. 
Equation  (20)  is  inappropriate  for  fitting  sueh  behavior,  and  in  Appendix  II  we  develop  a 
method  to  determine  a  lateney  time  to.  Any  population  information  that  oeeurred  before 
to,  in  the  groups  of  trajectories  that  displayed  latent  effeets,  was  excluded  from  the  fitting 
process  of  Eq.  (20). 

C.  Results  and  Discussion 

We  have  run  trajectories  on  all  three  PESs  to  compare  the  HO2  IVR  dynamics,  the 
isomerization  dynamics  of  H-0-0’  ^  0-0 ’-H,  and  the  unimoleeular  dissoeiation 
dynamies  of  HO2  ^11  +  02.  The  results  of  eaeh  of  these  are  now  discussed  in  order. 

l.IVR 

Normal  mode  kinetie  energies  as  fraetions  of  the  total  are  plotted  as  a  function  of  time 
in  Eigs.  7  and  8  for  total  energies  of  35.5  kcahmol  and  47.3  kcal/mol,  respectively;  the 
red,  magenta,  and  green  eurves  correspond,  respeetively,  to  the  results  of  the  initially 
excited  0-0  stretching,  0-0-H  bending,  and  0-H  stretching  modes.  These  two  different 
energies  were  seleeted  to  determine  the  effect  of  isomerization  on  the  relaxation  proeess. 
Isomerization  is  not  energetieally  accessible  on  any  of  the  PESs  at  35.5  kcal/mol  but  is  at 
47.3  kcal/mol.  The  results  in  the  columns  from  left  to  right  are  for  the  IMES,  XXZLG, 
and  DMBE-IV  PESs. 

The  results  for  35.5  kcal/mol  are  quite  similar  for  the  IMES  and  XXZEG  PESs  down 
to  fine  detail.  The  relaxation  timescale  is  well  beyond  the  duration  of  the  integration  time 
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of  the  trajectories  (10  ps)  regardless  of  the  excited  mode.  The  relative  ordering  of  the 
modes  from  fastest  to  slowest  is:  OH*,  OOH*  and  then  00*.  The  ordering  is  similar 
with  the  results  from  the  DMBE  IV;  however,  the  relaxation  of  trajectories  on  this 
surface  completed  within  1  ps  to  2  ps  (note  the  expanded  x-axis  timescale  for  the  DMBE- 
IV  PES  results  in  Eig.  7). 

The  results  in  Eig.  7  illustrate  some  of  the  limitations  of  our  mode-specific  trajectory 
studies.  The  ZPE  fraction  of  the  total  energy  of  35.5  kcal/mol  is  ~0.05  for  the  00  and 
OOH  modes  and  =0.15  for  the  OH  mode  (see  Table  111).  Consequently,  the  top  (00*) 
and  middle  (OOH*)  row  of  panels  in  Pig.  7  should  have  the  excited  mode  energy  fraction 
starting  at  =0.8  while  the  lower  row  (OH*)  should  start  at  =0.9.  This  expectation  of  the 
fractional  partitioning  of  energy  in  the  normal  modes  is  only  true  for  the  00*  energy. 
Por  the  other  cases,  the  initial  fraction  is  noticeably  lower,  and  for  OOH*  on  the 
DMBE-IV  PES,  it  would  appear  that  two  modes  are  excited.  On  all  three  PESs,  when 
either  OH*  or  OOH*  are  initially  excited,  there  is  a  rapid  energy  redistribution  from  the 
selected  initial  distribution  to  what  is  seen  in  Pig.  7  on  a  timescale  less  than  100  fs,  too 
fast  to  be  resolved  in  the  figure.  This  redistribution  is  undoubtedly  due  to  the  inherent 
difficulties  of  decomposing  the  total  energy  into  contributions  from  specific  modes  and 
the  use  of  normal  modes  for  assigning  the  initial  conditions.  The  kinetics  model  does  not 
include  such  processes,  and  no  attempt  is  made  to  initialize  the  model  to  the  true  initial 
conditions.  Rather,  as  discussed  in  the  Appendix  II,  a  least-squares  procedure  is  used  to 
get  effective  initial  conditions.  The  resulting  fits  are  the  black  lines  in  Pig.  7,  with 
optimized  values  of  ki,  k2,  and  appearing  c  in  Table  IV.  There  is  later  discussion  of 
expanding  the  two  exponentials  as  linear  functions  to  characterize  the  fitting  process. 
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These  linear  expansions  are  deseribed  by  two  slopes  sA  and  sT  that  also  appear  in  this 
Table  IV  . 

The  results  in  Fig.  7  show  that  the  kineties  model  fits  the  trajeetory  results  very  well. 
The  only  elear  but  minor  fitting  defieieney  oeeurs  for  the  IMLS  and  XXZLG  PESs  for 
OH*  at  early  times  (<  1  ps).  In  this  case,  the  trajectory  results  and  the  model  display  a 
fast  and  a  slow  decay  for  each  mode  of  energy.  However,  in  the  region  of  slow  decay  the 
00  and  OOH  energy  fractions  maintain  a  near  constant  absolute  difference  with  time,  a 
result  that  cannot  be  represented  by  the  common  slow  rate  constant  of  the  model.  In 
effect,  the  trajectory  results  exhibit  more  decay  constants  than  the  model  can  provide, 
implying  the  physical  principles  upon  which  the  model  is  based  are  insufficient  to 
describe  the  results  in  detail. 
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Figure  7  Ensemble-averaged  normal  mode  energies  as  functions  of  time  for  initial 
conditions  of  mode-specific  excitation  (see  text  for  details)  and  a  total  energy  of  =35.5 
kcaEmol.  From  the  top,  rows  are  results  for  the  excited  mode  being  00*,  OOH*,  and 
OH*.  From  the  left,  columns  are  results  for  the  IMFS,  XXZFG,  and  DMBE-IV  PESs.  In 
all  panels,  red,  magenta,  and  green  curves  represent  the  00,  OOH,  and  OH  modes, 
respectively,  while  the  black  lines  represent  the  kinetics  model  fits. 
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Table  IV.  IVR  rate  constants  (ps'*)  for  specific  mode  excitations. ‘‘ 


E  =35.5  kcal/mol 

E=47.3  kcal/mol 

IMLS 

XXZLG 

DMBE  IV 

IMES 

XXZEG 

DMBE  IV 

00* 

c  0.8(6) 

-0.1(3) 

1.66(9) 

2.2(9) 

2.1(3) 

4.2(1) 

ki 

0.003(5) 

0.017(4) 

0.72(5) 

0.05(3) 

0.07(1) 

0.32(1) 

ki 

0.014() 

0.037(4) 

2.18(5) 

0.16(3) 

0.22(1) 

1.37(1) 

sA 

-0.0076 

-0.012 

sr 

0.0024 

0.0045 

00* 

c  -2.0(7) 

1.2(4) 

-2.4(9) 

0.33(5) 

0.42(2) 

0.50(3) 

ki 

0.042(8) 

0.015(8) 

3.0(1) 

0.051(9) 

0.040(8) 

0.0(3) 

ki 

0.092(8) 

0.045(8) 

9.0(1) 

0.258(9) 

0.394(8) 

10.6(3) 

sA 

-0.021 

-0.018 

sr 

0.00075 

0.0045 

00* 

c  1.3(4) 

1.2(6) 

0.53(2) 

0.57(2) 

0.58(2) 

0.61(3) 

ki 

0.02(1) 

0.02(1) 

3.4(6) 

0.50(3) 

0.32(2) 

7.0(2) 

ki 

0.054(8) 

0.049(3) 

94.0(2) 

5.4(1) 

3.21(4) 

62.0(2) 

sA 

-0.032 

-0.025 

sr 

-0.00087 

0.0096 

^  Number  in  parentheses  is  the  uncertainty  in  the  last  digit. 
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Figure  8  is  the  analogue  of  Fig.  7  but  for  47.3  keal/mol,  whieh  is  above  the  barriers 
for  isomerization  for  the  three  PESs.  Table  IV  contains  the  optimized  values  of  the 
kinetics  model  parameters.  The  IVR  above  the  isomerization  threshold  is  similar  in  fine 
detail  for  the  IMLS  and  XXZLG  PESs.  The  relaxation  timescale  for  all  excited  modes  is 
comparable  to  or  less  than  the  duration  of  the  trajectories  (10  ps).  In  contrast,  on  the 
DMBE-IV  PES  IVR  is  complete  by  ~1  ps  or  less  no  matter  what  mode  is  excited  (notice 
the  smaller  x-axis  scale  for  the  DMBE-IV  PES  results  in  Eig.  8).  The  excited  OH*  and 
OOH*  relax  faster  with  00*  noticeably  slower  for  all  three  PESs,  as  is  the  case  for  the 
lower-lenergy  results  shown  in  Eig.  7.  As  in  Eig.  7  and  for  the  same  reasons,  there  are  fast 
redistribution  effects.  The  -12.3  keal/mol  increase  into  each  excited  mode  should  make 
the  initial  energy  fraction  higher  in  Eig.  8  than  in  Eig.  7.  In  almost  all  cases  the  apparent 
initial  fraction  is  lower  due  to  short  time  (<  100  fs)  energy  redistribution  not  resolvable  in 
the  panels. 
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Figure  8  The  same  as  Fig.  3  exeept  for  a  total  energy  of  =47.3  keal/mol. 

Although  the  results  in  Fig.  8  show  more  varied  relaxations  including  an  overshoot  of 
the  asymptotic  equipartition  of  energy,  the  kinetics  model  fits  the  trajectory  results  at 
nearly  comparable  quality  to  that  in  Fig.  7  for  the  lower  energy.  The  clearest  deficiency 
concerns  OOH*  relaxation  on  both  the  IMLS  and  XXZLG  PESs  where  the  decay  at  short 
times  is  underestimated.  In  this  case,  E2,  the  OFl  energy,  is  essentially  at  its  asymptotic 
value  of  E/3  except  at  times  too  early  to  be  resolved  on  the  plot.  Energy  conservation 
then  forces  Ei  and  £5  to  be  essentially  mirror  images  of  each  other  about  the  E2  at  all 
times.  The  fits  and  trajectory  results  display  this  symmetry  but  the  kinetics  model  is  too 
simplistic  to  represent  the  higher  frequency  or  transient  decay  components  of  the 
trajectories.  A  second  deficiency  is  OOH*  for  DMBE  where  E2  is  poorly  represented  at 
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very  short  times  (<200  fs).  The  temporal  range  of  the  panel  is  only  1  ps  but  the  fit  is 
applied  over  the  entire  10  ps  range.  Upon  close  examination  of  the  panel,  the  fit  has  two 
asymptotic  values,  one  for  E2  and  the  other  for  Ei  and  E3.  Through  most  of  the  range  of 
the  panel  and  beyond  to  the  10  ps  maximum  range  of  the  fit,  the  trajectory  results  for 
OOH*  on  the  DMBE-IV  PES  exhibit  an  E2  energy  fraction  persistently  lower  than  the 
energy  fraction  for  the  other  two  modes.  This  behavior  forces  the  optimum  fit  to  display 
two  asymptotes  at  the  price  of  a  poor  representation  of  E2  at  very  early  times.  (How  the 
model  can  realize  two  asymptotes  is  discussed  in  Appendix  II  material.)  In  this  case,  the 
trajectory  approach  to  equilibrium  is  too  complicated  for  the  kinetics  model. 

It  is  convenient  to  discuss  the  IMES  and  XXZEG  results  in  Table  IV  first,  then  the 
DMBE-IV  results.  A  leading  feature  of  the  IMES  and  XXZEG  results  is  that  for  35.5 
kcahmol,  all  the  values  for  kj  and  k2  are  smaller  than  0.1  ps"'  whereas  for  47.3  kcakmol  at 
least  k2  is  substantially  greater  than  0.1  ps'\  This  means  that  for  all  or  at  least  most  of  the 
10  ps  temporal  range  of  the  IVR  data  at  35.5  kcal/mol,  the  exponentials  in  the  fitting 
function  of  Eq.  (S.4),  found  in  Appendix  II,  can  be  reasonably  well  approximated  by  a 
first-order  expansion  that  is  linear  in  time.  The  fits  in  Eig.  7  clearly  show  this 
approximately  linear  behavior.  As  discussed  in  Appendix  II,  in  such  circumstances  the 
decay  is  best  represented  by  two  slopes,  5a  and  sr,  which  only  for  IMES  and  XXZEG  at 
35.5  kcal/mol  are  listed  in  Table  IV.  Eor  both  IMES  and  XXZEG,  the  values  of  5a  and  sr 
show  overall  that  IVR  becomes  more  rapid  in  the  order  00*,  OOH*,  and  OH*. 
However,  the  IMES  and  XXZEG  values  are  not  that  similar.  Eor  the  larger  slope  s^,  the 
IMES  value  varies  from  ~60%  lower  to  ~30%  higher  than  the  XXZEG  value.  The  lesser 
slope  Sr,  the  IMES  value  varies  from  =15%  to  90%  of  the  XXZEG  value.  The  IMES  and 
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XXZLG  results  at  35  kcal/mol  look  quite  similar  in  Fig.  7,  but  that  is  somewhat  due  to 
the  faet  that  only  a  fraction  of  the  decay  takes  place  in  10  ps.  Had  the  trajectory  results 
been  carried  out  to  longer  times,  the  differences  would  likely  be  more  manifest. 

The  IMLS  and  XXZLG  results  for  47.3  kcal/mol  in  Table  IV  show  relaxation  rate 
constants  that  are  all  considerably  larger  than  the  rate  constants  at  the  lower  energy. 
Exponential,  as  opposed  to  linear,  decay  is  much  more  evident  in  Fig.  8.  As  for  35 
kcahmol,  for  both  IMLS  and  XXZLG  the  IVR  decay  rate  constants  increase  in  the  order 
00*,  OOH*,  and  OH*.  The  OH*  rate  constants  are  an  order  of  magnitude  larger  than 
the  rate  constants  of  OOH*,  a  much  larger  increase  than  at  35  kcal/mol  and  an  indication 
of  an  accelerating  OH*  decay.  For  the  large  decay  rate  the  IMLS  value  is  from  =65% 
lower  to  =65%  higher  than  the  XXZLG  value.  For  the  smaller  decay  rate  kj,  the  IMLS 
value  is  from  =70%  lower  to  =60%  higher  than  the  XXZLG  value.  Together  with  similar 
results  at  35  kcal/mol,  the  results  indicate  that  IMLS  and  XXZLG  do  have  somewhat 
different  IVR  decay  characteristics.  Over  this  energy  range,  IMLS  IVR  occurs  more 
slowly  for  00*  and  more  rapidly  for  OH*  relative  to  XXZLG  IVR  decay.  For  OOH* 
over  this  energy  range,  the  IMLS  and  XXZLG  exchange  order  with  respect  to  decay  rate 
constants. 

Unlike  IMLS  and  XXZLG,  a  10  ps  timescale  is  sufficient  at  both  energies  to  capture 
all  of  the  IVR  processes  on  DMBE-IV.  The  ki  and  k2  values  time  constants  obtained 
using  the  DMBE-IV  PES  in  Table  IV  are  at  least  40  times  larger  at  35  kcal/mol  and  at 
least  4.5  times  larger  at  47.3  kcal/mol  than  the  corresponding  rate  constants  computed 
using  the  IMLS  or  XXZLG  PESs.  This  implies  a  significant  qualitative  difference  in  the 
dynamics  and  associated  IVR.  For  example,  the  DMBE-IV  results  for  OH*  at  both 
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energies  below  and  above  isomerization  has  the  larger  k2  time  constant  exceeds  50  ps'\ 
implying  relaxation  over  within  a  10  fs  to  20  fs,  or  one  to  two  OH  vibrational  periods. 
Such  processes  are  not  resolvable  in  the  IMLS  and  XXZLG  panels  of  Figs.  7  and  8  but  in 
the  expanded  timescale  of  the  panels  for  DMBE,  they  are  just  visible  in  the  earliest  times 
on  the  plots.  As  discussed  earlier,  on  this  timescale  limitations  in  the  initial  condition 
preparation  and  the  approximations  used  to  monitor  mode-specific  energies  lead  to 
noticeable  apparent  redistributions  of  energy.  For  IMLS  and  XXZLG,  the  fit  is 
insensitive  to  such  processes  because  slow  energy  transfer  processes  that  extend  over  the 
whole  temporal  range  outweigh  such  brief  transients.  However,  for  DMBL-IV,  these 
transients  are  more  comparable  to  the  usual  energy  transfer  processes  and  the  fits  at  least 
of  OH*  reflect  that. 

Consistent  with  IMLS  and  XXZLG,  the  DMBL-IV  decay  rate  is  slowest  for  00*  and 
fastest  for  OH*.  In  contrast  to  IMLS  and  XXZLG,  the  DMBL-IV  decay  rate  constants 
have  a  much  milder  dependence  on  energy  and,  in  the  case  of  00*,  the  decay  rate 
constants  actually  decrease  with  energy.  Careful  examination  of  the  00*  results  for  the 
DMBL-IV  PES  in  Ligs.  7  and  8  show  the  trajectories  predict  a  slightly  slower  relaxation 
at  the  higher  total  energy.  The  DMBL-IV  PES  has  the  highest  isomerization  barrier,  ~1.8 
kcahmol  higher  than  IMLS.  This  might  imply  that  the  phase  space  volume  associated 
with  the  isomerization  region  of  the  PES  is  less  accessible  for  the  DMBL-IV  PES  than 
for  the  IMLS  or  XXZLG  PES,  resulting  in  a  corresponding  decrease  to  already  rapid  IVR 
for  the  bend  and  stretch  modes  most  coupled  at  the  isomerization  saddle  point. 

To  better  understand  the  effect  of  isomerization  on  IVR,^°^  the  energy  fraction  of 
different  modes  is  plotted  as  a  function  of  time  in  Lig.  9  for  the  time-dependent  ensemble 
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of  trajectories  that  have  not  yet  isomerized.  As  time  progresses,  isomerization  reduces  the 
number  of  trajectories  in  this  ensemble  causing  the  sampling  noise  in  the  ensemble 
averaged  energy  fraction  to  increase.  For  this  reason,  the  energy  fraction  was  not 
followed  beyond  the  time  where  90%  of  the  trajectories  had  isomerized.  The  comparison 
of  Fig.  9,  which  excludes  isomerization,  to  Fig.  8,  which  includes  isomerization,  shows 
qualitative  differences. 
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Figure  9  The  same  as  Fig.  8  only  for  the  populations  of  those  trajectories  at  time  t  that 
have  not  yet  isomerized.  The  time  history  of  each  unisomerized  population  stops  at  either 
10  ps,  or  when  only  10%  of  the  population  remains. 

Before  discussing  the  differences,  note  that  the  OOH*  and  OH*  DMBE-IV  results  in 
panels  (h)  and  (i)  are  almost  identical  in  Figs.  8  and  9.  As  will  be  shown  in  Sec.  C.2  of 
the  results  and  discussion,  for  these  two  cases  the  isomerization  rate  is  significantly 
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slower  than  the  IVR  rate.  The  IVR  is  essentially  eomplete  before  isomerization  begins. 
In  all  other  eases  in  Figs.  8  and  9,  the  IVR  is  eomparable  to  or  slower  than  isomerization. 
In  these  other  eases,  Fig.  9  shows  that  isomerization  seleetively  removes  from  the 
ensemble  those  trajeetories  that  spend  most  of  their  time  with  more  OOH  bend  energy 
and  leaves  behind  trajeetories  that  have  larger  00  and  OH  stretch  energy.  This  is  most 
clearly  seen  for  all  three  PESs  for  00*.  The  results  in  Fig.  9  suggest  that  the  asymptotic 
equipartition  of  energy  is  not  reached  in  the  ensemble  of  unisomerized  radicals  because 
the  isomerization  removes  those  with  sufficient  energy  in  the  OOH  and  OH  modes  to 
isomerize  leaving  a  population  with  “excess”  energy  in  the  00  mode.  As  shown  in 
Sec.  B.2,  the  main  difference  between  DMBE-IV  and  IMES  or  XXZEG  in  the  top  row  of 
Eig.  9  is  that  the  isomerization  and  IVR  rate  constants  are  comparable  for  only 
DMBE-IV,  for  which  some  of  the  IVR  is  complete  before  isomerization  becomes 
significant.  In  the  cases  of  OOH*  and  OH*  excitations  on  the  IMES  or  XXZEG  PES,  the 
removal  of  isomerized  species  from  the  ensemble  results  in  a  much  shorter  timescale 
shown  in  Eig.  9,  which  is  an  indication  of  how  dominant  isomerization  is  over  IVR.  The 
results  for  OOH*  clearly  illustrate  that  isomerization  is  much  faster  than  IVR;  note  that 
the  timescale  of  Eig.  9  is  about  one-tenth  that  of  Eig.  8.  An  important  difference  in  the 
two  PESs  is  that  they  differ  by  ~2  kcal/mol  in  the  isomerization  barrier. 

The  results  in  Eig.  9  show  how  isomerization  acts  to  create  a  new  energy  distribution 
in  the  as  yet  unisomerized  population.  Under  normal  circumstances,  that  new  energy 
distribution  cannot  be  accessed  because  eventual  re-isomerization  would  further  mix  all 
populations.  However,  conceptually  if  the  isomerized  HO2  could  be  prevented  from  re- 
isomerizing,  then  the  filtering  by  isomerization  could  have  chemical  consequences.  Eor 
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example,  if  two  different  isotopes  of  oxygen  were  present  in  HO2  and  if  some  isotopieally 
seleetive  photoehemical  process  could  prevent  a  re-isomerization,  then  the  unisomerized 
population  could  carry  out  chemistry  with  its  filtered  energy  distribution.  This  kind  of 
effect  would  be  present  for  any  molecule  where  isomerization  is  faster  than  IVR. 

2.  Isomerization 

Since  the  trajectories  used  for  the  IVR  calculations  in  Figs.  8  and  9  isomerized  on  all 
three  PESs,  these  same  trajectories  can  be  analyzed  for  mode-specific  isomerization.  The 
normalized  populations  of  unisomerized  trajectories  as  a  function  of  time  are  plotted  in 
Fig.  10  along  with  the  fit  according  to  Eq.  (20),  whose  rate  constants  are  listed  in  Table 
V.  In  the  figure,  red  (00*),  magenta  (OOH*),  and  green  (OH*)  curves  are  compared  to 
the  black  dashed  curves  of  the  fit. 

Figure  10  displays  what  we  have  called  a  latency  time  in  Sec.  B.2  For  all  three  PESs, 
00*  has  an  initial  accelerating  rate  while  for  IMFS  and  XXZFG,  OH*  has  a  shelf  visible 
in  the  plot  inserts.  A  prominent  reason  for  such  features  is  that  the  initial  position  of  the 
hydrogen  atom  for  00*  and  OH*  are  well  removed  from  the  perpendicular  bisecting 
plane  of  the  0-0  bond.  Even  with  facile  energy  transfer  into  the  bending  mode,  the 
hydrogen  atom  will  take  time  to  reach  the  isomerization  region.  Also  as  discussed 
concerning  IVR,  imperfect  initial  conditions  can  affect  initial  energy  transfer.  Taken 
together,  these  effects  can  give  rise  to  latency  in  the  isomerization  time.  As  described  in 
Supplemental  Material,  the  population  up  to  the  latency  time  listed  in  Table  V  is  not 
included  in  the  fit  of  Eq.  (20).  This  results  in  a  single  exponential  fit,  i.e.,  5  =  0  in  Eq. 
(20),  for  all  three  excitations  for  DMBE-IV  and  for  00*  for  IMES  and  XXZEG. 
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Figure  10  For  the  IMLS  (a),  XXZLG  (b),  and  DMBE-IV  (c)  PESs,  the  population  decay 
as  a  function  of  time  due  to  isomerization  for  initial  conditions  of  mode-specific 
excitation  at  a  total  energy  of  -47.3  kcahmol  (see  text  for  details).  In  all  panels,  red, 
magenta,  and  green  curves  represent  the  00,  OOH,  and  OH  modes  respectively  while  the 
black  lines  represent  the  kinetics  model  fits.  Inserts  in  each  panel  magnify  the  results  at 
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early  times. 

A  prominent  reason  for  sueh  features  is  that  the  initial  position  of  the  hydrogen  atom 
for  00*  and  OH*  is  well  removed  from  the  perpendicular  bisecting  plane  of  the  0-0 
bond.  Even  with  facile  energy  transfer  into  the  bending  mode,  the  hydrogen  atom  will 
take  time  to  reach  the  isomerization  region.  Also  as  discussed  concerning  IVR,  imperfect 
initial  conditions  can  affect  initial  energy  transfer.  Taken  together,  these  effects  can  give 
rise  to  latency  in  the  isomerization  time.  As  described  in  Appendix  II,  the  population  up 
to  the  latency  time  listed  in  Table  V  is  not  included  in  the  fit  of  Eq.  (20).  This  results  in  a 
single  exponential  fit,  i.e.,  5  =  0  in  Eq.  (20),  for  all  three  excitations  for  DMBE-IV  and 


for  00*  for  IMES  and  XXZEG. 
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Table  V.  Isomerization  rate  eonstants  (ps'^),  prefactors,  and  latency  times  (ps)  under  different  initial  conditions  as  a 
function  of  energy  for  the  three  different  PESs. 


Energy 

(kcal/mol) 

Initial 

Conditions 

IMES 

A  ki 

B 

ki 

to 

XXZEG 

A  ki 

B 

ki 

to 

DMBE-IV 
A  ki 

to 

44.0 

EMS 

0.56 

0.04 

0.40 

0.71 

0.49 

0.13 

0.48 

1.31 

1.00 

0.22 

46.0 

EMS 

0.46 

0.16 

0.52 

1.77 

0.40 

0.19 

0.54 

1.93 

1.00 

0.53 

47.3 

00* 

0.99 

0.13 

1.53 

0.99 

0.18 

1.24 

1.23 

0.68 

0.38 

OOH* 

0.08 

0.23 

0.87 

3.40 

0.14 

1.04 

0.84 

5.01 

0.92 

0.77 

OH* 

0.16 

0.31 

0.99 

1.63 

0.17 

0.11 

0.30 

1.07 

1.98 

0.28 

0.94 

0.82 

48.0 

EMS 

0.43 

0.21 

0.55 

2.53 

0.37 

0.32 

0.59 

2.73 

1.00 

0.90 

50.0 

EMS 

0.33 

0.28 

0.62 

2.84 

0.27 

0.36 

0.69 

2.94 

0.96 

1.44 

52.0 

EMS 

0.30 

0.36 

0.66 

3.49 

0.30 

0.58 

0.68 

4.14 

0.97 

1.78 
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For  IMLS  and  XXZLG,  Fig.  10  clearly  shows  the  rate  constant  for  isomerization  is 
quite  dependent  on  which  mode  was  initially  excited.  In  particular,  the  00* 
isomerization  rate  is  on  the  order  of  10  times  slower  than  the  rate  when  either  of  the  two 
other  modes  are  exeited.  In  Table  IV,  the  relevant  IVR  rate  k2  is  eomparable  to  the 
isomerization  rate  in  Table  V.  Since  0-0  vibration  is  not  the  isomerization  reaetion 
coordinate,  IVR  is  a  necessary  first  step  for  isomerization  to  oeeur.  The  substantial 
latency  times  and  the  slow  isomerization  rate  indicate  that  IVR  is  the  bottleneck  for 
isomerization  for  00*.  As  would  be  expeeted,  OOH*  has  the  fastest  isomerization  rate 
and  no  latency  time  because  bending  motion  is  a  dominant  component  of  the 
isomerization  reaction  path.  The  dominant  isomerization  rate  is  an  order  of  magnitude 
faster  than  the  IVR  rate  constants  in  Table  IV,  meaning  that  isomerization  is  largely  over 
before  substantial  amounts  of  energy  have  leaked  by  IVR  into  modes  less  well  connected 
to  the  isomerization  reaetion  path.  The  slower  rate  in  the  bi-exponential  fit  to  OOH* 
isomerization  is  somewhat  eomparable  to  the  IVR  rate  constants  in  Table  IV,  suggesting 
that  this  ~10%  tail  to  the  declining  population  is  a  consequence  of  IVR  processes  re¬ 
assembling  energy  into  the  bend.  The  OH*  isomerization  rate  eonstants  fall  in  between 
those  of  00*  and  OOH*.  From  Table  IV,  the  fastest  IVR  rate  constants  are  two  or  three 
times  faster  than  the  dominant  isomerization  rate  constants  in  Table  V  and  comparable  to 
the  dominant  isomerization  rate  constants  for  OOH*.  Although  faster  than  isomerization, 
there  is  a  time  lag  between  converting  OH  stretching  motion  into  OOH  bending  motion 
which  leads  to  the  small  (~100  fs)  latency  times  listed  in  Table  V.  IVR  distributes  energy 
to  OOH,  promoting  isomerization,  but  also  to  00,  which  is  inefficiently  connected  to 
isomerization.  Like  OOH*,  the  ~10%  tail  to  the  population  deelines  at  a  slow  rate  that  is 
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comparable  to  most  of  the  IVR  rate  eonstants  in  Table  IV,  suggesting  an  IVR-type 
bottleneck  to  re-assembling  initially  misdirected  energy  baek  into  the  bend. 

In  contrast  to  IMLS  and  XXZLG,  the  results  in  Fig.  10c  and  Table  V  clearly  show  the 
isomerization  rate  obtained  using  DMBE-IV  is  independent  of  which  mode  is  excited. 
The  main  differenee  between  modes  seen  in  Fig.  10c  is  not  the  rate  but  the  latency  time 
for  00*.  Contrasting  Tables  IV  and  V,  the  00*  faster  IVR  rate  is  twice  that  of 
isomerization,  leading  to  a  latency  time  smaller  by  a  factor  of  three  than  those  of  the 
other  two  PESs.  The  00*  isomerization  rate  is  the  smaller  than  that  from  either  OOH* 
or  OH*  excitations  by  ~15%,  perhaps  due  to  some  lingering  effects  of  IVR  processes. 
The  OOH*  and  OH*  faster  IVR  rate  constants  are  at  least  an  order  of  magnitude  faster 
than  the  isomerization  rate  constants,  meaning  that  for  most  of  isomerization,  memory  of 
the  initial  excitation  is  nearly  lost.  There  are  no  lateney  times  and  the  isomerization  rate 
constants  are  within  10%  of  eaeh  other. 

The  isomerization  deeay  of  an  initially  random  population  is  plotted  in  Eig.  1 1  for 
total  energies  44,  46,  48,  50,  and  52  kcal/mol  on  eaeh  PES. 
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Figure  11  For  the  IMLS  (a),  XXZLG  (b),  and  DMBE-IV  (c)  PESs,  the  population  decay 
as  a  function  of  time  due  to  isomerization  for  initial  conditions  selected  using  EMS  (J=0) 
sampling  and  for  total  energies  in  kcal/mol  of  44  (red),  46  (magenta),  48  (green),  50 
(brown),  and  52  (cyan).  The  black  lines  represent  the  kinetics  model  fits. 
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This  range  of  energies  was  seleeted  so  the  radieal  would  be  above  the  respeetive 
thresholds  of  isomerization  for  all  of  the  PESs,  but  below  the  elassieal  thresholds  for 
dissoeiation  (see  Table  I).  In  Figs.  11(a)  and  11(b),  the  deeay  eurves  for  IMLS  and 
XXZLG  are  fit  with  a  bi-exponential  function  (shown  by  the  dotted  curves  in  the  figure) 
to  produce  the  fast  and  slow  rate  constants  in  Table  V  for  the  grid  of  total  energies.  In 
contrast,  the  decay  curves,  computed  using  the  DMBE-IV  PES  results,  are  well  fit  by  a 
single  exponential;  the  values  of  the  parameters  from  the  fits  are  given  in  Table  V.  The 
bi-exponential  rate  constants  of  isomerization  for  statistical  initial  conditions  can  be 
interpreted  by  using  the  graphical  representation  in  Fig.  12  of  the  rate  constants  in  Table 
V  plotted  as  a  function  of  the  total  energy  above  the  isomerization  barrier  on  each  PES 
(see  Table  I).  The  results  are  color  coded  by  PES:  red,  IMES;  blue,  XXZEG;  magenta, 
DMBE-IV.  The  two  solid  straight  lines  at  the  bottom  of  the  plot  are  the  linear  least 
squares  fits  to  the  slower  EMS  isomerization  rate  constants  for  IMES  and  XXZEG.  The 
three  dotted  straight  lines  in  the  figure  are  the  linear  least-squares  fits  to  the  faster  IMES 
and  XXZEG  EMS  isomerization  rate  constants  and  the  single  DMBE-IV  statistical 
isomerization  rate  constants.  The  open  symbols  (squares,  triangles,  diamonds)  are  in 
order  the  00*,  OOH*,  and  OH*  isomerization  rate  constants  for  DMBE-IV  or  IMES  and 
XXZEG  (only  the  rate  with  the  largest  prefactor).  (The  OOH*  XXZEG  faster 
isomerization  rate  of  5.01  ps'^  the  largest  rate  in  Table  IV,  is  off  the  scale  of  the  plot.) 
The  total  energy  in  excess  of  the  isomerization  barrier  is  used  for  the  x  axis  in  this  plot  to 
make  a  side-by-side  comparison  of  the  rate  constants  clearer. 
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Figure  12  Rate  constants  for  isomerization  from  Table  V  plotted  as  functions  of  the  total 
energy  above  the  isomerization  barrier  (Viso).  See  text  for  details 

For  IMLS  and  XXZLG,  among  the  lowest  isomerization  rate  constants  are  those  for 
00*.  These  rate  constants  are  somewhat  lower  than  the  straight-line  fits  to  the  slower 
EMS  isomerization  rate  constant.  The  likely  interpretation  is  that  the  slow  rate  on  either 
IMLS  or  XXZLG  is  due  to  a  bottleneck  in  relaxing,  from  the  initial  statistical 
distribution,  those  trajectories  that  have  a  high  degree  of  0-0  vibrational  excitation.  The 
calculated  rate  is  a  reflection  of  both  isomerization  and  energy  transfer  from  the  0-0 
stretch  into  the  two  other  degrees  of  freedom  that  are  coupled  to  the  isomerization  saddle 
point.  In  contrast,  among  the  highest  IMLS  or  XXZLG  rate  constants  are  those  for  either 
OOH*  or  OH*.  At  the  same  energy,  the  OH*  rate  constants  are  comparable  to  the  faster 
EMS  rate  constants  while  the  OOH*  rate  constants  are  approximately  double  the  faster 
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EMS  rate  constants  (note  the  OOH*  XXZLG  rate  is  off  the  scale  of  the  plot). 
Consequently,  the  faster  EMS  rate  constants  on  IMES  and  XXZEG  have  no  significant 
energy  transfer  bottleneck  from  either  the  0-H  stretch  or  the  0-0-H  bend  modes.  These 
faster  rate  constants  are  a  closer  measure  of  isomerization  uncomplicated  by  energy 
transfer  issues.  Eor  DMBE-IV,  as  discussed  above,  the  00*,  OOH*,  and  OH* 
isomerization  rate  constants  are  nearly  the  same  because  IVR  is  fast  enough  to  scramble 
any  mode  specific  initial  conditions.  Since  EMS  initial  conditions  are  statistically 
random,  it  is  not  surprising  that  the  least-squares  fits  to  the  EMS  isomerization  rate 
constants  essentially  pass  through  the  mode-specific  rate  constants.  The  00* 
isomerization  rate  is  the  lowest  below  the  linear  fit  to  the  statistical  rate  constants  and  in 
this  case,  as  discussed  above,  the  IVR  rate  is  only  double  the  isomerization  rate. 
Interestingly,  despite  the  different  topologies  seen  in  the  three  PESs  in  Eig.  5,  the 
component  of  the  IMES,  XXZEG,  and  DMBE-IV  isomerization  rate  constants  least 
contaminated  by  IVR  effects  all  have  a  similar  size,  within  a  factor  of  two,  as  indicated 
by  the  dotted  lines  in  Eig.  12.  This  occurs  only  if  the  barrier  height  is  subtracted  from  the 
total  energy,  as  in  Eig.  12,  indicating  the  typical  dominance  of  barrier  height  over  saddle 
point  shape  is  determining  the  gross  features  of  the  rate  constant. 

3,  Unimolecular  dissociation 

The  mode  specific  behavior  found  in  IVR  and  isomerization  motivates  a  similar 
investigation  for  the  unimolecular  dissociation  of  HO2  ^  H  +  O2.  As  discussed  earlier, 
only  the  XXZEG  and  DMBE-IV  PESs  are  fully  specified  at  this  limit  and  can  be  studied. 
Marques  and  Varandas  carried  out  exactly  such  a  study  on  the  DMBE-IV  PES,  portions 
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of  which  we  repeat  with  more  trajeetories  and  extend  to  XXZLG.  As  discussed  in  See. 
C.l  we  used  fractional  vibrational  quantum  numbers  on  the  XXZLG  PES  in  order  to 
preserve  the  same  total  energies  above  Dg.  The  approximate  total  energy  listed  in  Table 
III  is  the  average  energy  for  all  six  mode-speeific  runs.  While  the  energy  in  any  exeited 
mode  is  always  below  the  dissociation  energy,  the  energy  of  the  OH*  mode  is  especially 
close  to  the  dissoeiation  energy  being  only  =0.1  keal/mol  and  =0.2  kcal/mol  below  the 
dissoeiation  limits  for  the  DMBE-IV  and  XXZEG  PESs,  respeetively.  In  addition  to 
mode-specific  initial  conditions,  Marques  and  Varandas^°  also  studied  statistieal  initial 
conditions  at  one  dissociative  energy,  whieh  we  repeat  and  extend  to  the  XXZEG  PES  at 
a  shifted  total  energy  to  preserve  the  excess  energy  above  Dg.  The  resulting  unimolecular 
deeay  eurves  are  shown  in  Eig.  13  and  the  rate  constants  derived  from  them  in  Table  VI. 
The  DMBE-IV  results  are  consistent^®^  with  the  results  of  Marques  and  Varandas. 
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Table  VI  For  HO2  ^  H  +  O2,  unimolecular  dissociation  rate  constants  (ps'^),  prefactors,  and  latency  times  (ps)  under  different  initial 
conditions  as  a  function  of  the  total  energy  for  the  two  different  PESs.  The  energy  in  parenthesis  is  the  amount  of  energy  in  excess  of 
the  De  in  Table  I. 


Initial 

conditions 

Energy 

(kcahmol) 

DMBE-IV 

A  ki 

B 

k2  to 

Energy 

(kcal/mol) 

XXZEG 

A  ki 

B 

k2 

to 

00* 

58.79(3.95) 

1.18 

0.69 

0.29 

57.63(3.95) 

1.18 

0.47 

0.41 

OOH* 

58.74(3.90) 

1.00 

0.76 

57.57(3.89) 

0.79" 

0.49" 

0.25" 

1.04" 

0.10 

OH 

58.26(3.41) 

0.86 

0.67 

57.09(3.41) 

0.47 

0.45 

0.41 

12.5 

EMS 

59.43(4.59) 

1.00 

0.89 

58.27(4.59) 

0.87 

0.62 

0.10 

3.32 

The  least  squares  solution  shows  an  rms  fitting  error  that  is  a  weak  function  of  correlated  small  changes  in  ki  and  large  changes  in 
k2.  Many  more  trajectories  would  be  needed  to  determine  k2  with  a  precision  comparable  to  other  entries  in  the  table. 
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Figure  13  For  the  XXZLG  (a)  and  DMBE-IV  (b)  PESs,  the  population  deeay  as  a 
funetion  of  time  due  to  unimoleeular  dissoeiation  to  H  +  O2  for  initial  eonditions  of 
mode-speoifie  excitation  [00*  (red),  OOH*  (magenta),  OFl*  (green)]  or  for  initial 
conditions  selected  using  EMS  (J=0)  sampling  (brown).  The  black  lines  represent  the 
kinetics  model  fits.  Inserts  in  each  panel  magnify  the  results  at  early  times.  See  text  for 
details. 

Three  points  about  the  results  in  Eig.  13  and  Table  VI  are  worth  noting.  Eirst,  in  Eig. 
13  on  both  PESs  there  are  significant  early  time  behaviors  that  cause  the  mode-specific 
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population  decays  to  generally  be  distinct  from  each  other.  The  early  time  decay  for  each 
excitation  mode  will  now  be  considered  in  the  order  00*,  OH*,  and  OOH*. 

On  both  PESs,  00*  has  the  most  pronounced  latency  effects  as  detected  by  the 
procedure  outlined  in  Appendix  I.  Comparisons  of  Fig.  13  and  Table  VI  for  dissociation 
with  Fig.  4  and  Table  V  for  isomerization  shows  that  the  latency  effects  in  00* 
dissociation  are  similar  in  form  but  shorter  than  those  for  00*  isomerization  on  both 
PESs.  Since  insufficiently  fast  IVR  rate  constants  are  a  likely  source  of  this  latency  for 
isomerization,  can  IVR  00*  rate  constants  at  dissociation  energies  be  estimated?  For 
XXZFG,  the  linear  extrapolation  of  the  00*  IVR  rate  constants  in  Table  IV  produce  a 
rate  at  the  00*  dissociation  energy  of  Table  VI  that  is  =0.4  ps'\  less  than  but  comparable 
to  the  00*  dissociation  rate,  a  situation  analogous  to  that  in  isomerization.  This 
extrapolation  is  bolstered  by  the  fact  that,  as  discussed  above,  the  slower  EMS 
isomerization  rate  is  thought  to  track  the  00*  IVR  rate.  As  seen  in  Fig.  12,  this  slow  rate 
is  reasonably  well  approximated  as  a  linear  function  of  energy.  Extension  of  this  linear 
function  to  dissociation  energies  requires  an  extrapolation  over  half  the  energy  range 
required  with  extrapolation  using  Table  IV.  The  end  result  is  a  slow  rate  component  of 
isomerization  that  is  higher  than  but  comparable  to  the  dissociation  rate.  Thus  the  latency 
in  00*  dissociation  likely  has  the  same  cause  as  the  latency  in  isomerization  with  00*. 
The  latency  in  00*  dissociation  is  shorter  than  that  for  isomerization  because  both  IVR 
and  dissociation  rate  constants  are  higher  at  the  higher  energy.  For  DMBE-IV, 
extrapolations  from  Table  IV  are  dubious  because,  as  commented  above,  the  IVR  rate 
constants  in  Table  IV  decrease  with  energy,  a  result  most  unlikely  to  persist  at  higher 
energies.  Furthermore,  there  is  no  particular  reason  to  correlate  00*  IVR  rate  constants 
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on  the  DMBE  IV  with  EMS  isomerization  rate  constants.  However,  unlike  the  case  for 


XXZEG,  the  results  in  Tables  V  and  VI  show  that  the  00*  isomerization  rate  at  the 
higher  energy  at  which  IVR  was  calculated  is  almost  identical  to  the  00*  dissociation 
rate.  Thus  if  the  IVR  rate  constants  were  to  just  flatten  out  with  energy  or  modestly 
increase,  the  relationship  of  IVR  rate  constants  to  dissociation  rate  constants  would  be 
comparable  to  that  for  isomerization.  Since  IVR  is  faster  for  DMBE-IV  than  XXZEG,  for 
both  isomerization  and  dissociation  the  latency  time  is  shorter  for  DMBE-IV  than  for 
XXZEG.  Since  this  scenario  does  not  involve  a  significant  increase  in  any  of  the  rate 
constants  from  that  for  isomerization,  the  latency  time  at  dissociation  is  closer  to  that  for 
isomerization  than  is  the  case  for  XXZEG. 

Eor  OH*,  there  is  no  latency  time  on  either  PES.  However,  as  the  inserts  in  Eig.  13 
show,  ~10%  of  the  population  has  the  last  inner  turning  point  in  negative  time.  In  other 
words,  10%  of  the  trajectories  are  ballistic,  that  is  they  are  created  with  initial  conditions 
that  directly  lead  to  dissociation. As  seen  in  Tables  I  and  III,  the  OH*  excitation 
energy  is  only  =0.2  kcahmol  below  the  dissociation  energy  and  trajectories  need  to  gain 
only  a  few  tenths  of  a  kcal/mol  from  the  other  two  modes  to  dissociate.  The  OH*  results 
on  the  DMBE-IV  PES  suggest  that  trajectories  that  have  one  inner  turning  point  are  much 
more  likely  to  lose  significant  energy  out  of  the  OH  stretch  and  begin  a  statistical 
repopulation  of  the  phase  space  of  the  dissociative  coordinate.  However,  the  OH*  results 
on  the  XXZEG  PES  suggest  that  a  few  turning  points  do  not  significantly  decrease  the 
energy  in  the  OH  stretch,  allowing  trajectories  during  each  expansion  phase  of  the  OH 
stretch  to  acquire  from  the  two  other  modes  the  small  amount  of  energy  they  need  to 
dissociate.  The  difference  between  the  two  PESs  is  undoubtedly  a  reflection  of  the  order 
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of  magnitude  higher  OH*  IVR  rate  for  DMBE-IV  relative  to  XXZLG  seen  in  Table  IV 
for  lower  total  energies.  This  difference  is  reflected  in  the  fit  to  the  two  different 
dissociation  decay  curves.  For  DMBE-IV,  the  fitting  error  is  not  sensitive  to  the  very 
early  time  behavior,  so  although  the  local  fitting  error  is  an  order  of  magnitude  larger 
than  at  any  other  time,  the  ~5  fs  duration  of  this  very  large  local  error  is  overwhelmed  by 
minimizing  the  global  error  over  the  remaining  10  ps  of  decay.  For  related  reasons,  this 
initial  short  period  is  not  recognized  as  a  latency  time  by  the  procedure  outlined  in  Sec. 
B.2  However,  for  XXZEG  because  the  ballistic-like  mechanism  persists  much  longer,  bi¬ 
exponential  behavior  results  with  the  very  rapid  process  having  a  rate  constant  of  ~13  ps"' 
(see  Table  VI).  Since  the  vibrational  period  of  the  0-H  stretch  is  ~10  fs,  the  inverse  of 
this  fast  rate  ~15  fs  implies  that  for  half  a  dozen  vibrational  periods  the  OH*  retains 
sufficient  energy  for  the  radical  to  undergo  ballistic  dissociation  before  enough  energy  is 

1  AT 

lost  out  of  the  OH  stretch  to  require  the  statistical  reassembly  of  sufficient  energy  to 
dissociate. 

For  OOH*,  there  is  no  latency  on  the  DMBE-IV  PES.  This  is  as  expected  because 
any  reasonable  extrapolation  of  the  OOH*  IVR  results  in  Table  IV  to  dissociation 
energies  would  have  OOH*  IVR  more  than  an  order  of  magnitude  larger  than  the  OOH* 
dissociation  rate  in  Table  VI.  However,  for  XXZEG,  doing  the  same  would  have  IVR 
rate  constants  somewhat  higher  than  but  comparable  to  the  OOH*  dissociation  rate 
constants.  This  would  suggest,  based  on  the  arguments  in  the  00*  case,  a  noticeable 
latency  time  for  dissociation.  There  is  in  fact  a  latency  time  of  only  60  fs,  which  is  too 
short  to  influence  the  OOH*  dissociation  rate.  The  reason  for  such  a  small  latency  was 
discussed  in  connection  with  Figs.  7  and  8.  These  IVR  figures  show  plots  of  mode 
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specific  energies  as  a  function  of  time,  but  the  initial  distributions  in  the  figures  show 
much  less  OOH  energy  and  much  more  OH  energy.  The  explanation  is  a  massive 
redistribution  of  the  initial  mode-specific  energy  in  tens  of  femtoseconds  that  cannot  be 
distinguished  within  the  time  resolution  of  the  figures.  This  is  particularly  prominent  at 
the  higher  energy  (see  Fig.  8)  where  the  OH  mode  appears  to  start  with  the  asymptotic 
energy. 

This  short  time  redistribution  is  due  at  least  in  part  to  the  breakdown  in  the  normal 
mode  approximation  in  the  initial  conditions  preparation.  However,  the  normal  mode 
approximation  is  more  valid  for  the  0-0-H  bending  potential  than  it  is  for  the  stretches. 
For  the  stretches,  compression  leads  to  a  rapid  increase  in  repulsive  energy  not  well 
described  by  the  harmonic  potential  of  normal  modes.  In  the  bending  potential,  both 
compression  and  expansion  lead  to  a  barrier  through  a  symmetry  plane  (see  Fig.  5)  on  all 
three  PESs.  In  HO2,  the  isomerization  barrier  is  lower  than,  but  comparable  to,  the  barrier 
through  the  collinear  geometry.  Consequently,  the  harmonic  potential  is  a  somewhat 
better  approximation  of  the  anharmonic  potential  than  it  is  for  the  stretches. 

Beyond  the  breakdown  of  the  normal  mode  approximation,  the  ballistic  mechanism 
discussed  above  for  OH*  suggests  a  second  very  rapid  energy  transfer  mechanism.  A 
careful  examination  of  the  early  time  OOH*  populations  shown  in  the  inserts  in  Fig.  10 
reveals  a  behavior  for  OOH*  isomerization  similar  to  that  seen  in  the  inserts  in  Fig.  13 
for  OH*  dissociation.  A  ballistic  interpretation  of  these  isomerization  population  features 
would  have  been  clearer  if  a  turning  point  definition  in  terms  of  the  last  0-0-H  bend 
turning  point  had  been  used,  because  that  would  have  led  to  negative  lifetimes  for 
rigorously  ballistic  trajectories  as  is  the  case  for  OH*  dissociation.  Isomerization  might 
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very  well  accelerate  energy  transfer  out  of  OOH*.  The  massive  redistribution  of  energy 
in  the  OOH*  fVR  case  occurs  on  the  same  timescale  as  ballistic  isomerization,  namely  a 
few  oscillations  of  the  bending  mode.  This  redistribution  is  much  more  pronounced  in 
Fig.  8  for  energies  above  dissociation  than  in  Fig.  7  for  energies  below  dissociation.  This 
redistribution  is  also  somewhat  more  pronounced  for  OOH*  than  for  OH*.  Finally,  as 
seen  in  Figs.  7  and  8,  the  primary  beneficiary  of  loss  of  energy  out  of  OOH*  is  the  OH 
mode  which  is  the  reaction  coordinate  for  dissociation.  Energy  flow  between  bending  and 
stretching  modes  can  be  especially  facile;  see,  for  example.  Refs.  24(a)  and  108.  Thus 
the  breakdown  in  the  normal  mode  approximation  and  ballistic  isomerization  could  both 
contribute  to  very  early  IVR  and  unexpected  short  dissociation  latency  times  for  OOH* 
on  XXZLG.  Furthermore,  OOH*  dissociation  on  XXZLG  obeys  bi-exponential  behavior 
like  that  for  OH*  dissociation  (see  Table  VI).  The  larger  secondary  rate  constant  for 
OOH*,  although  much  smaller  than  the  larger  OH*  rate  constant  and  with  half  the 
amplitude,  may  be  due  to  ballistic  trajectories.  In  the  case  of  OOH,*  a  probable  source  of 
the  ballistic  dissociative  trajectories  may  well  be  ballistic  isomerization  which,  as  long  as 
it  lasts,  initially  pumps  energy  into  the  OH  reaction  coordinate  to  accelerate 
dissociation. 

The  second  conclusion  that  can  be  drawn  from  Fig.  13  and  Table  IV  is  that  both  PESs 
display,  at  best,  modest  mode-specific  effects  in  unimolecular  dissociation.  While  Eig. 
13  has  distinct  curves  depending  on  the  mode  excited,  especially  for  XXZLG,  those 
differences  primarily  reflect  the  latency  and  ballistic  effects  at  early  times  discussed 
extensively  above.  As  shown  in  Table  VI,  the  dominant  rate  constants  from  the  mode 
specific  conditions  are  quite  similar  to  each  other  with  only  an  8%  (5%)  maximum 


73 


90 


deviation  from  the  average  for  DMBE-IV  (XXZLG).  If  plotted  as  funetions  of  energy, 
these  mode-speeifie  unimoleeular  dissoeiation  rate  eonstants  on  both  surfaees  are 
reasonably  eonsistent  with  the  EMS  dissoeiation  rate  eonstants  at  somewhat  higher 
energies.  In  the  ease  of  isomerization,  mode-speeifie  effeets  for  XXZEG,  espeeially  for 
00*,  are  prominent  and  eomparison  with  the  IVR  rate  eonstants  ealeulated  at  the  same 
energy  suggests  eompetition  between  IVR  and  isomerization.  Eor  OH*  and  OOH* 
dissoeiation  on  both  surfaees,  IVR  is  generally  fast  enough  to  not  be  a  substantial 
bottleneek  to  dissoeiation.  Eor  00*,  IVR  at  short  times  is  slow  enough  to  eause  lateney 
effeets  on  both  PESs.  However,  for  DMBE-IV,  Eig.  12  would  suggest  that  the 
isomerization  rate  at  the  00*  dissoeiation  energy  will  be  on  the  order  of  5  times  higher 
than  the  00*  dissoeiation  rate  in  Table  VI.  Eor  XXZEG,  as  diseussed  previously,  the 
slower  eomponent  of  the  EMS  isomerization  rate  is  assoeiated  with  the  00* 
isomerization  rate  and  at  the  00*  dissoeiation  energy  it  would  be  at  least  double  the  00* 
dissoeiation  rate  in  Table  VI.  To  the  degree  that  isomerization  drives  energy  transfer,  the 
IVR  rate  at  early  times  that  influenees  lateney  may  not  be  the  IVR  rate  operational  over 
the  bulk  of  the  unimoleeular  dissoeiation  after  the  lateney  time. 

The  third  and  final  eonelusion  is  that,  as  ean  be  seen  from  Table  VI,  DMBE-IV 
produees  unimoleeular  dissoeiation  rate  eonstants  that  are  approximately  50%  larger  than 
does  the  XXZEG  at  the  same  exeess  energy  over  the  dissoeiation  energy.  There  have 
been  four  previous  studies  ’  ’  ’  that  eompared  reaetions  on  these  two  PESs  but  all  of 
them  were  quantum  dynamies  studies  of  reaetions  (la)  and  (lb),  not  the  unimoleeular 
reaetion  (le)  eonsidered  here.  The  more  relevant  studies^^’^^  with  H  +  O2  as  reaetants 
predieted  a  lower  reaetion  probability  for  XXZEG  relative  to  DMBE-IV.  The  two  other 


74 


91 


studies^^’^'^  with  O  +  OH  as  reactants  found  that  only  at  intermediate  energies  was 
XXZLG  more  reactive.  The  results  in  Table  VI  are  thus  consistent  with  the  general 
observation  that  XXZLG  is  less  reactive  than  DMBE-IV.  The  quantum  dynamics 
studies  also  make  clear  that  the  reaction  dynamics  of  HO2  is  dominated  by  quantum 
mechanical  resonances  unrepresented  in  classical  trajectories  that  can  result  in  large 
discrepancies  in  classical  microcanonical  and  even  canonical  rate  constants  relative  to 
their  quantum  analogs. 

One  last  point  concerning  Table  VI  is  the  relationship  between  two  bi-exponential 
cases  in  the  table:  OH*  and  EMS  dissociation  on  XXZEG.  (The  XXZEG  OOH*  case  is 
also  bi-exponential  but  not  to  a  degree  sufficient  to  influence  the  following  discussion.) 
The  OH*  case  is  described  above  as  a  combination  of  an  early  time  ballistic-type 
dissociation  process  and  the  usual  statistical  dissociation  process.  Most  likely  the  EMS 
behavior  is  bi-exponential  because  it  contains  a  faint  hint  of  the  OH*  ballistic  process.  If 
one  assumes  that  the  EMS  distributed  energy  is  approximated  by  an  equal  mix  of  00*, 
OOH*,  and  OH*  initial  conditions,  then  the  quite  similar  slow  rate  constants  of  00*, 
OOH*,  and  OH*  would  be  reflected  in  a  similar  EMS  rate  that  has  increased  due  to  the 
higher  total  energy,  and  more  particularly  the  relatively  higher  excess  energy,  for  the 
EMS  case.  However  the  fast  OH*  rate  constants  should  also  be  reflected  in  a  fast  EMS 
rate,  with  only  1/3  weight  to  the  prefactor  if  we  assume  that  each  mode  contributes 
equally,  or  =0.13.  In  fact  the  amplitude  of  the  fast  EMS  rate  is  =0.10.  This  fast  rate, 
while  =5  times  faster  than  the  slow  rate,  is  also  about  3  times  slower  than  the  OH*  fast 
rate.  These  differences  are  probably  due  to  both  absorbing  the  latency  features  of  the 
00*  dissociation  and  the  severity  of  approximating  EMS  by  equal  portions  of  00*, 
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OOH*,  and  OH*.  While  this  approximation  is  indeed  severe,  the  EMS  case  produces 
some  trajectories  with  last  inner  turning  points  prior  to  dissociation  occurring  at  negative 
times.  Even  though  the  energy  is  statistically  distributed,  some  trajectories  do  have  OH 
energies  comparable  to  the  dissociation  energy  as  the  above  interpretation  suggests. 

D.  Conclusions 

Three  PESs  for  HO2  have  been  comparatively  examined  for  their  characterization 
of  equilibrium,  isomerization,  and  dissociation  energies,  structures,  and  quantum 
variational  vibrational  frequencies,  phase  space  structures  using  surfaces-of-section,  and 
for  IVR,  isomerization,  and  dissociation  trajectory  dynamics.  Two  of  the  PESs,  DMBE- 
and  have  been  fully  described  in  the  literature  while  the  third 

surface,  IMES,  has  been  only  introduced  in  the  literature^^  and  is  more  fully  described 
here.  The  two  most  recent  PESs,  XXZEG  and  IMES,  are  based  on  higher  quality  and 
more  systematically  distributed  ab  initio  electronic  structure  calculations  than  the 
DMBE-IV,  a  PES  for  HO2  that  has  been  used  in  many  studies.  The  XXZEG  and  IMES 
differ  in  the  manner  in  which  the  electronic  structure  calculations  are  fit  and  a  slight 
improvement  in  the  quality  of  the  ab  initio  method  in  going  from  XXZEG  to  IMES. 
These  two  PESs  are  not  without  limitations  with  IMES  not  being  globally  characterized 
up  to  the  lowest  dissociation  channel  and  XXZEG  displaying  small  violations  of 
permutation  symmetry  and  some  fitting  imperfections  due  to  data  sparsity. 

The  major  conclusions  of  this  study  are: 

•  The  IMES  and  XXZEG  are  quite  similar  and  of  high  quality.  In  comparison  to 
experiment,  the  rms  error  of  IMES  and  XXZEG  vibrational  frequencies  is  ~I0  cm"' 
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while  the  dissoeiation  energy  is  too  low  by  -0.2  kcal/mol  for  IMLS  and  =0.8  kcal/mol 
for  XXZLG.  Surfaces-of-seetion  plots  as  a  funetion  of  energy  show  consistent  and 
comparable  phase  space  structures  for  trajectories  calculated  on  both  the  IMLS  and 
XXZLG.  Irregular  phase  space  structure  is  first  apparent  at  energies  near  the 
isomerization  barrier.  The  structured  nature  of  phase  space  gives  rise  to  slow  and 
roughly  comparable  IVR  processes  on  both  PESs.  The  relaxation  of  excited  00  stretch 
is  especially  slow.  Slow  IVR  in  turn  leads  to  prominent  mode-specific  effects  in 
isomerization  (on  both  PESs)  and  dissociation  (only  accessible  on  XXZEG).  Since 
IVR  is  not  fast  enough  to  redistribute  the  initial  vibrational  energy,  even  for  the  initial 
conditions  selected  by  using  EMS  isomerization,  decay  curves  display  double 
exponential  behavior  with  the  slower  rate  correlated  with  slow  IVR  rate  constants.  The 
likely  source  of  this  slow  IVR  is  energy  initially  localized  in  the  00  stretch.  The 
longer  it  takes  for  a  trajectory,  with  sufficient  energy,  to  isomerize;  the  more  probable  it 
becomes  that  energy  has  been  localized  in  the  00  mode.  As  the  00  mode  becomes 
progressively  enriched  with  energy  the  chance  for  isomerization  decrease.  In  other 
words,  isomerization  acts  as  a  filter  that  over  time  that  can  leave  behind  a  population 
highly  excited  in  the  00  stretching  mode. 

•  In  detail  XXZEG  and  IMES  do  display  somewhat  different  dynamics.  The  heights  of 
the  isomerization  saddle  points  on  the  two  PESs  differ  by  =2.0  kcal/mol.  (This  is 
primarily  due  to  sparsity  in  the  XXZEG  ab  initio  database  around  the  symmetry  plane 
of  the  surface  and  not  differences  in  the  underlying  electronic  structure  calculations). 
Results  sensitive  to  isomerization  have  an  energy  dependence  shifted  approximately  by 
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that  amount.  The  IVR  rate  constants  computed  on  the  two  PESs  differ  by  a  factor  of 
two  for  some  initial  conditions. 

•  The  DMBE-IV  PES  displays  many  different  characteristics  relative  to  IMES  and 
XXZEG.  It  does  not  compare  well  to  experiment  with  ~100  cm"'  rms  errors  in 
vibrational  frequencies  and  -0.7  kcal/mol  overestimation  of  the  dissociation  energy. 
The  surfaces-of-section  plots  show  a  more  irregular  phase  space  structure  starting  at 
energies  that  are  half  the  isomerization  barrier  height.  The  IVR  decay  rate  constants  are 
nearly  an  order  of  magnitude  faster  than  those  for  IMES  and  XXZEG.  As  a 
consequence,  mode-specific  effects  on  isomerization  and  unimolecular  dissociation  are 
much  less  prominent  and  single  exponential  decay  indicative  of  statistical  behavior  was 
found. 

•  The  isomerization  and  unimolecular  dissociation  rate  constants  were  determined  with  a 
kinetics  model  that  identified  early  time  latency  effects  by  comparing  the  fitting  error  to 
the  sampling  error  due  to  the  finite  number  of  trajectories  in  each  simulation.  This 
approach  works  reasonably  well.  In  a  related  fashion,  the  IVR  decay  on  all  three  PESs 
was  analyzed  with  a  kinetics  model  based  on  an  asymptotic  limit  of  equipartition  of 
energy  and  on  mode-to-mode  energy  transfer  rate  constants  that  are  proportional  to  the 
amount  of  energy  in  the  initial  mode.  This  model  also  works  well  and  represents  the 
IVR  features  on  all  three  PESs  at  energies  above  and  below  isomerization.  However, 
the  derived  mode-to-mode  rate  constants  were  sometimes  negative  and  therefore  not 
physically  meaningful,  reducing  the  model  to  a  phenomenological  one.  This  may  be 
due  to  the  effects  of  resonances*^*^'^^’^^^  on  the  energy  transfer  that  are  not  accounted  for 
in  the  model. 
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•  Dissociation  rate  constants  could  only  be  determined  on  XXZLG  and  DMBE-IV.  On 
both  surfaces,  latency  times  due  to  slow  energy  transfer  and  fast  ballistic  processes 
result  in  population  decays  that  are  distinctly  different.  However,  after  these  early  time 
effects,  the  dominant  dissociation  rate  constants  are  at  best  only  modestly  sensitive  to 
initial  mode-specific  excitations.  As  is  true  for  bimolecular  reactions,  XXZLG  has  a 
lower  dissociation  rate  constant  than  DMBE-IV  at  comparable  energies  above  the 
dissociation  limit. 

On  balance  the  IMLS  and  XXZLG  are  the  PESs  of  choice  with  IMLS  slightly 
superior  because  of  better  energetics  and  a  more  reliable  description  of  isomerization. 
However,  the  IMLS  at  present  cannot  be  used  for  dynamics  studies  at  energies  above 
dissociation  without  incorporating  additional  ab  initio  information. 

E.  Directions  of  Future  Research 

There  are  three  directions  of  research  that  could  be  conducted  on  the  dynamics  and 
reactivity  of  HO2.  The  first  direction  is  investigating  the  quantum  mechanical  effect  of 
tunneling  which  cannot  be  accounted  for  due  to  the  classical  nature  of  these  studies.  The 
effect  of  proton  tunneling  would  be  manifest  by  vibrational  line  splitting  and  by  an 
enhanced  rate  of  isomerization.  There  have  been  similar  studies  done  with  the  hydrogen 
peroxyl  anion.  The  PESs  of  the  anion  and  the  radical  share  some  similar  topological 
characteristics,  such  as  the  symmetric  hydrogen  wells  and  a  barrier  to  isomerization.  On 
the  anion  PES  the  classical  barrier  to  isomerization  is  ~17  kcal/mol,  which  is  a  little 
under  half  the  barrier  height  of  the  radical.  At  the  isomerization  saddle  point  the  0-0 
bond  length  is  elongated  relative  to  the  bond  distance  at  the  minima  as  it  is  with  the 
radical.  The  higher  barrier  to  isomerization  on  the  radical  could  result  in  tunneling  being 
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less  important  at  lower  energies.  Previous  studies  of  proton  transfer  studies  with  the 
anion  found  mode  specific  effects  with  excitation  of  the  0-0-H  mode  contributing  to  the 
largest  splitting  factors  and  there  was  mode  specific  enhancement  from  excitation  of  the 
0-0  mode.  It  would  also  be  of  interest  to  see  if  the  lower  isomerization  barrier  with  the 
anion  would  lead  to  enhanced  IVR  compared  to  the  IVR  found  with  the  radical. 

The  second  direction  of  future  research  with  HO2  is  related  to  the  reactivity  of  the 
radical.  The  study  conducted  so  far  only  looked  at  the  dissociation  of  the  radical  in  a 
very  narrow  energy  range,  but  in  this  energy  range  modest  mode  specific  effects  were 
found  and  the  lifetime  of  the  radical  was  bi-exponential.  A  general  study  of  the 
dissociation  reactions,  across  a  larger  energy  range,  could  give  more  insight  into  the 
nonstatistical  dynamics  of  the  radical.  Such  trajectory  studies  would  be  complemented 
by  comparison  to  predictions  from  statistical  rate  calculations.  Monte  Carlo  variational 
transition  state  theory  (MCVTST),  ^^°given  the  barrierless  nature  of  the  H-0  bond  fission, 
would  be  useful  for  comparison  since  the  study  can  be  conducted  on  the  same  potential 
energy  surface  as  the  trajectory  study,  which  would  facilitate  a  direct  comparison.  A  rate 
in  MCVTST  theory  is  calculated  using^"^*^^^’^^'^*^^ 


KE)  =  \ 


\drdmE)-E)d{q,c-qMq„c 


(21) 


\dr5{H{T)-E) 

V  / 

where  F  is  the  phase  space  of  the  system,  H(r)  is  the  Hamiltonian  with  no  motion  of  the 
center  of  mass,  qRc  is  the  reaction  coordinate  which  is  a  function  of  the  other  coordinates 


q,  qc  is  the  value  required  for  reaction  or  a  dividing  surface  and  ^j-^is  the  flux  through  qc 
MCVTST  rate  constants  have  previously  been  interpreted  as  the  upper  limit  to  rate 
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constants  calculated  by  trajectories.  This  means  that  the  inequality  k(E)McvTST  <  k(E)traj 
would  be  observed  if  the  trajectories  did  not  sample  all  of  the  energetieally  accessible 
phase  space.  Sueh  a  comparison  would  be  useful  in  further  understanding  the  dynamics 
ofH02. 

We  investigated  the  energy  transfer  between  the  vibrational  modes  of  HO2  on 
three  different  PESs  and  found  evidence  of  different  timeseales  of  relaxation  whieh  was 
shown  to  be  dependent  on  the  potential  energy  surface  used. ' ' '  However,  in  these  studies 
the  radical  was  not  rotating.  The  final  proposed  direction  of  future  research  would  focus 
on  understanding  the  coupling  between  rotation  and  vibration.  Frederiek  et  al.  used 
trajeetories  to  study  OCS  and  SO2  and  found  that  large  seale  energy  transfer  between 
rotation  and  vibration  (R-V)  readily  occurred.  It  is  probable  that  centrifugal  effects  from 
rotation  motion  could  influence  the  rate  constants  of  isomerization  and  dissociation  in 
HO2. 
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IV.  Dynamics  of  Hydrogen  Peroxyl  Radical  in  a  Dense  Gas 
Environment 

A.  Introduction 

Understanding  intermolecular  energy  transfer  is  an  important  step  in  the  meehanistic 
understanding  of  H2/O2  eombustion.  We  are  interested  in  the  collisional  deactivation  of 
HO2,  because  the  transfer  of  energy  into  and  out  of  the  internal  motions  of  the  radical  will 
ultimately  control  whether  the  radical  stabilizes  or  reacts.  If  the  energized  radical 
decomposes,  it  propagates  the  free  radical  chain;  if  it  is  stabilized  by  collisions,  the  free 
radical  chain  terminates.  Kinetics  modeling  of  HO2  in  combustion  environments  requires 
mapping  out  a  large  number  of  interconnected  reaction  pathways  and  the  associated 
reaction  and  energy  transfer  rates,  but  possible  inaccuracies  in  the  time  constants  used  in 
parameterizing  combustion  models  can  result  in  disagreements  between  theoretical 
predictions  and  experimental  results.  Completely  resolving  all  of  the  discrepancies 
between  theory  and  experiment  is  a  monumental  task,  so  in  hopes  of  contributing  useful 
information  about  the  rich  chemistry  of  this  radical,  we  continue  with  a  “from  the  bottom 
up”  approach  on  HO2.  In  chapter  III,  which  has  been  published,^''  we  focused  on  the 
intramolecular  dynamics  of  the  radical  and  found  complex  dynamics  and  multiple 
timescales  for  energy  transfer  amongst  the  vibrational  modes  of  the  radical,  which  affect 
the  unimolecular  dissociation  of  the  radical.  In  a  similar  spirit  we  next  focus  on  studying 
the  deactivation  of  a  vibrationally  excited  HO2  in  a  dense-gas  environment. 

We  have  chosen  to  simulate  vibrationally  excited  HO2  in  an  argon  bath.  Argon  was 
selected  as  the  solvent  bath  to  avoid  any  possible  phase  change  of  the  medium,  since  the 
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Ar  bath  was  simulated  at  the  temperatures  and  pressures  above  the  experimental  eritieal 
point/ The  timeseales  for  relaxation  in  a  dense  gas  or  liquid  environment  are 
eommonly  diseussed  in  the  eontext  of  two  physieal  proeesses:  Ti  for  vibrationally 
deaetivating  eollisions  and  T2  for  elastie  collisions  that  dephase  the  motion  of  the 
solute/'^  The  isolated  binary  collision  approximation  (IBC)  is  commonly  used  to 
calculate  the  timescale  and  behavior  of  T/  as  a  function  of  density.  Litovitz^^^  was  one  of 
the  first  to  apply  a  gas  phase  binary  collision  model  in  describing  relaxation  of  an  excited 
species  in  a  liquid.  The  gist  of  this  theory  is  that  relaxation  of  a  vibrationally  excited 
species  in  a  fluid  can  be  described  as  a  series  of  isolated  collision  events  that  are 
uncorrelated  from  one  another.  The  effect  of  density  on  the  energy  transfer  process  can 
be  described  by  a  simple  multiplication  of  an  energy  transfer  probability  with  the 
collision  frequency.  The  energy  transfer  probability  is  assumed  to  be  independent  of  the 
internal  energy  of  the  solute.  The  collision  frequency  is  calculated  for  the  particular 
density  or  pressure  of  the  system;  at  the  simplest  approximation  the  collision  frequency  is 
based  on  hard  sphere  model. 

Early  debates  concerning  the  validity  of  IBC  theory  centered  on  neglecting  the 
possible  effects  of  many  body  collisions.  Herzfeld  pointed  out  that  infrequent,  high 
energy  collisions  would  contribute  the  most  to  energy  transfer  into  and  out  of  the 
vibrational  degrees  of  freedom,  so  solute  relaxation  consists  of  a  large  number  of  elastic 
collisions  punctuated  with  occasional  inelastic  collisions.  Davis  and  Oppenheim^^^ 
revised  the  definition  of  IBC  by  using  a  pairwise  radial  distribution  function  (RDF)  to 
implicitly  include  many  body  effects  in  the  collision  definition.  Their  definition  of 
collision  frequency  is 
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V,  =  AnR^pgiR) 


/  \l/2 

kT 

^np) 


(22) 


where  47rR^  defines  a  spherical  shell  at  R,  p  is  the  system  density,  g(R)  is  the  RDF 
evaluated  at  R  and  [xT/itp)  is  the  average  thermal  velocity  representing  the  flux  of  the 
bath  gas  through  the  shell.  They  also  emphasized  that  the  only  important  collisions  are 
those  that  explore  high  energy  regions  of  the  potential,  along  the  same  lines  of  reasoning 
as  Herzfeld.  Dardi  and  Cukier  put  forth  an  interesting  criticism  that  the  selection  of  R 
in  Eq.  (22)  is  arbitrary  and  a  sensible  collision  radius  cannot  be  defined  for  both  low  and 
high  density  regions.  They  pointed  to  and  reanalyzed  the  experimental  work  of  Chatelet 
et  al.  on  the  relaxation  of  Hi  and  showed  that  the  relaxation  time  constants  at  higher 
densities  do  not  extrapolate  to  the  low-density  predictions. 


The  proposal  that  high  energy  collisions  are  the  dominant  source  of  relaxation  was 
explored  by  Ohmine,  using  simulations  of  excited  ethylene  in  both  Ar  and  HiO. 
Ohmine  found  for  Ar  that  high  energy  collisions,  which  explore  the  repulsive  wall  of  the 
potential,  dominate  the  relaxation,  but  multistep  vibration-to-rotation  (V-R)  energy 
transfer  followed  by  rotational  energy  transfer  into  the  bath  plays  a  role.  Relaxation  in 
water  was  found  to  involve  much  shorter  timescales,  which  Ohmine  attributed  to  the 
much  stronger  electrostatic  interactions  with  the  solute  and  the  rapid  response  of  the 
solvent  to  rearrangements  in  the  system.  Modifying  the  Ar-ethylene  parameters  to 
increase  the  magnitude  of  the  interaction  was  found  to  accelerate  the  relaxation. 

Schultz  et  al.  studied  the  relaxation  of  azulene  in  a  rare  gas  solvent  and  used  the 
Davis  and  Oppenheim  definition  of  collision  frequency  to  understand  the  weak  density 
dependence  of  the  relaxation.  They  found  that  the  implicit  assumption  of  spherical 


84 


101 


symmetry  in  this  collision  model  led  to  an  over-estimation  of  the  actual  density 
dependence  of  the  relaxation.  Using  MD  simulations  with  a  rigid  azulene  embedded  in  a 
dense  bath,  they  found  that  the  anisotropy  of  the  solute  and  the  intermolecular  potential 
led  to  a  density  independent  packing  above  and  below  the  azulene  rings.  They  illustrated 
this  anisotropy  by  calculating  RDFs  with  explicit  consideration  of  the  spatial  orientation 
of  the  solute.  These  RDFs  indicate  a  much  higher  local  density  in  the  region  above  and 
below  the  plane  of  the  molecule,  particularly  at  higher  densities.  They  suggest  that  the 
bath-gas  packing  could  inhibit  the  energy  transfer  out  of  the  azulene;  because,  the  lowest 
frequency  modes  in  the  molecule  involve  out-of-plane  motion  that  could  be  constricted 
by  the  packing  of  the  bath  gas  atoms. 

Schwarzer  et  studied  the  vibrational  relaxation  of  excited  azulene  in  several 
solvents  by  monitoring  the  decay  of  a  UV  absorption  band.  At  higher  pressures 
corresponding  to  liquid  densities  there  was  a  fast  initial  decay  of  the  signal,  which 
became  increasingly  exponential-like  at  longer  times.  They  found  a  nonlinear 
dependence  between  the  relaxation  times  of  azulene  and  solvent  density  in  Xe,  CO2,  and 
CiHg.  They  ruled  out  local  heat  conduction  because  aphysical  adjustments  to  the  model 
would  have  to  be  made  to  fit  the  data.  They  also  reasoned  that  IVR  was  unimportant 
since  the  deviations  from  linear  behavior  occurred  at  different  scaled  collision 
frequencies  for  the  various  bath  gases.  In  other  words,  they  reasoned  that  IVR  was  not  a 
limiting  factor  since  the  deviation  from  linear  behavior  did  not  occur  at  the  same  collision 
frequency.  They  proposed  a  simple  adsorption  model  of  azulene  forming  collision 
complexes  with  the  bath  gas  which  they  posited  would  decrease  the  collision  frequency 
and  consequently  inhibit  energy  transfer.  They  assumed  that  the  reduction  in  the 
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collision  frequency  would  be  linearly  dependent  on  the  fractional  coverage  of  the  solute 
by  the  bath  gas  atoms.  They  found  that  the  model  gave  agreement  at  lower  densities,  but 
failed  at  higher  densities.  They  postulated  that  a  ehange  in  the  average  energy  transferred 
per  collision  could  account  for  the  failure  of  the  model  at  higher  densities. 

Benzler  et  al.  did  a  similar  UV  absorption  study  to  measure  the  relaxation  of 
cycloheptatriene  and  found  a  nonlinear  density  dependence  on  the  relaxation  timescales 
as  the  solvent  density  inereased.  Along  a  similar  line  of  reasoning  as  used  by  Schwarzer 
et  they  concluded  that  local  heating  and  IVR  were  not  important  factors  in  the 

nonlinear  density  dependence  of  the  relaxation.  They  proposed  a  modified  version 
collision-eomplex  model  of  Sehwarzer  et  al.,  except  that  fractional  coverage  of  the  solute 
by  the  bath  gas  atoms  was  assumed  to  follow  exponential  dependence,  rather  than  a  linear 
one.  This  model  gave  acceptable  fits  to  the  experimental  data  and  they  noted  that 
saturation  of  the  speetral  shift  in  the  absorption  band  at  higher  pressure  could  be  related 
to  the  effect  of  bath  gas  atoms  clustering  around  the  solute,  but  additional  work  would 
be  needed  to  confirm  this  model. 

A  series  of  trajectory  studies  on  the  relaxation  of  exeited  azulene  embedded  in 
and  were  motivated  by  the  density  dependent  relaxation  time  constants 

observed  by  both  Sehwarzer  et  al.  and  Benzler  et  al.  The  first  study  by  Heildelbach 
et  a/.  consisted  of  simulations  at  two  different  pressures  and  two  different 
temperatures.  The  azulene  was  simulated  with  two  different  PES  models:  one  whieh 
reduced  the  number  of  intramoleeular  parameters  and  constrained  C-H  bonds  and  one 
which  represented  the  full  motion  of  the  intramolecular  potential.  The  CO2  moleeules 
were  also  represented  with  redueed  and  full  PESs,  which  were  used  with  the  respective 
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reduced  and  full  azulene  PESs.  The  authors  of  this  study  focused  on  investigating  gross 
relaxation  time  constants,  local  heating  effects,  collision  frequency  calculations  and  the 
local  structure  of  the  solvent.  The  second  trajectory  study was  directed  towards 
investigating  how  the  vibrational  modes  of  azulene  coupled  to  the  bath  and  how  the  time 
constants  of  energy  loss  occurred  on  a  mode -by-mode  basis. 

Heildelbach  et  al.  found  in  the  first  study  that  the  gross  relaxation  time  constants 
could  be  fit  to  a  single  exponential  for  both  the  low  and  high  pressure  simulations.  It 
should  be  noted  that  the  integration  time  for  the  lower  pressure  (250  atm)  is  ~3  times 
shorter  than  the  experimental  timescale  reported  by  either  Schultz  et  al}^^  or  Schwarzer 
et  The  high  pressure  results  were  simulated  for  ~l/8*  of  the  reported  experimental 
timescales  of  Schwarzer  et  al.  for  comparable  pressures.  The  short  timescales  of  these 
simulations  and  incomplete  relaxation  of  the  azulene  preclude  making  any  additional 
observations  about  the  complete  thermalization  process.  The  relaxation  was  heavily 
influenced  by  whether  full  or  reduced  PES  representations  of  the  azulene  and  CO2  were 
used.  No  evidence  of  local  heating  was  found  in  examining  the  kinetic  energy  of  the  bath 
atoms  in  close  proximity  to  the  solute.  The  structure  of  the  bath  around  the  solute  was 
found  to  be  anisotropic,  but  the  anisotropy  above  and  below  the  plane  of  the  azulene  was 
related  only  qualitatively  to  the  calculated  RDEs  and  the  explicit  orientation  of  azulene 
was  not  considered.  The  RDEs  calculated  by  Schultz  et  are  much  more  illustrative 
of  this  anisotropy.  The  formation  of  bath  gas-solute  collision  complexes,  which  were  an 
earlier  suggested  cause  of  nonlinear  density  dependence  of  relaxation  timescales,  was  not 
investigated 
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The  interaction  of  the  vibrational  modes  with  the  bath  was  the  motivation  for  the 
second  paper  with  the  focus  being  primarily  on  the  CO2  bath.  The  simulation  details 
were  similar  to  those  described  in  the  first  paper.  The  second  study  also  contained  details 
of  additional  simulations  of  azulene  embedded  in  Xe.  Heildelbach  et  al.  found  that 
energy  loss  from  azulene  in  low  pressure  CO2  occurred  primarily  through  the  low 
frequency  modes,  but  at  elevated  pressures  the  higher  frequency  modes  began  to  have  a 
larger  role  in  the  relaxation  process.  The  relaxation  in  Xe  occurred  primarily  through  the 
low  frequency  modes.  In  both  cases  the  low  frequency  motion  corresponding  to 
out-of-plane  modes  of  the  molecule  contributed  the  most  to  the  energy  transfer.  They 
also  found  that  the  translational  and  rotational  energy  of  the  molecule  relaxed  rather 
quickly  and  that  V-R  transfer  was  unimportant  in  the  relaxation.  They  also  studied  the 
IVR  of  isolated  azulene  using  both  simulated  frictional  forces  and  simple  non¬ 
equilibrium  distributions  of  energy  amongst  the  vibrational  degrees  of  freedom.  They 
found  that  the  relaxation  of  the  modes  in  the  isolated  molecule  happened  on  longer 
timescales  and  that  non-equilibrium  energy  distributions  could  occur.  They  proposed  that 
the  solvent  plays  a  role  in  facilitating  the  IVR,  which  ultimately  led  to  single  timescale 
for  energy  loss  to  the  bath. 

Paul  et  studied  the  relaxation  of  hexafiuorobenzene  embedded  in  a  bath  of  N2  at 
several  densities  using  classical  trajectories,  at  a  single  bath  temperature  of  298  K.  The 
initial  state  of  the  CeFe  was  vibrationally  hot.  The  translational  and  rotational  motion  of 
the  molecular  was  selected  to  be  thermal.  They  found  only  minimal  V-R  energy  transfer 
occurring  during  the  relaxation  process.  The  effect  of  the  intermolecular  potential  on 
relaxation  dynamics  was  examined  and  determined  to  have  minimal  effect.  Most  of  the 
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energy  lost  by  the  solute  was  transferred  to  the  translational  and  rotational  degrees  of 
freedom  in  the  bath.  The  relaxation  of  the  solute  was  well  fit  by  a  bi-exponential 
funetion  with  the  asymptotic  energy  being  treated  as  an  adjustable  parameter.  Based  on 
the  asymptotic  energy  of  CeFe  they  concluded  that  a  bi-exponential  was  best  suited  for 
lower  densities  and  that  additional  exponential  terms  would  be  needed  to  properly  model 
the  higher  density  results.  The  groups  of  96  trajectories  that  were  simulated  for  each 
density  were  decomposed  into  smaller  subgroups  for  analysis  of  the  number  of 
trajectories  needed  to  converged  results.  It  was  found  the  same  general  behavior  and 
fitting  results  could  be  obtained  with  a  sample  of  24  trajectories.  Interestingly,  it  appears 
that  Paul  et  al.  were  not  aware  that  the  nonlinear  density  dependence  of  the  relaxation 
time  constants  they  observed  is  quite  similar  to  the  trends  seen  in  the  results  of  Benzler  et 
al.  and  Schwarzer  et  al.  An  interesting  contrast  between  these  two  studies  is  that  the 
relaxation  of  CeFe  was  found  to  be  bi-exponential  instead  of  single  exponential. 

To  our  knowledge  there  are  only  two  studies  of  the  collisional  relaxation  of 
vibrationally  excited  HO2  and  both  are  based  on  bimolecular  trajectory  scattering.  The 
first  study,  by  Brown  and  Miller,  focused  on  the  deactivation  of  H02by  collisions  with 
He  as  a  function  of  bath  temperature,  total  molecular  energy,  and  total  angular 
momentum.  They  concluded  that  there  are  two  ways  that  HO2  is  vibrationally 
deactivated.  Low  temperature  collisions  (800  K)  with  a  rotationally  cold  radical  tended 
to  cause  an  internal  V-R  energy  transfer.  Hotter  collisions  (2,000  and  5,000  K)  with  the 
bath  gas  with  a  rotationally  excited  HO2  resulted  in  more  direct  loss  of  vibrational  energy 
to  translational  energy  of  the  radical. 
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The  second  study,  by  Varandas  and  Zhang, focused  on  the  deactivation  of  HO2  by 
collisions  with  O2.  The  relative  translational  temperature  of  O2  colliding  HO2  varied 
from  ~300  K  to  -4000  K,  and  the  rotational  and  vibrational  states  of  O2  were 
approximately  assigned  to  the  ground  state.  The  initial  vibrational  state  of  HO2  consisted 
of  three  vibrational  energies  (21,  36,  and  48  kcal/mol)  that  were  selected  by  quasi- 
classical  excitation.  The  rotational  temperature  of  HO2  was  -300  K.  They  found  the 
dominant  relaxation  mechanism  was  V-R  transfer,  across  all  of  the  impact  temperatures, 
and  that  collisions  generally  caused  a  V-R  energy  transfer.  A  mechanistic  observation 
shared  with  Brown  and  Miller.  They  also  found  that  the  shape  of  the  distribution  of 
energy  transferred  per  collision  is  only  weakly  dependent  on  the  initial  vibrational  energy 
of  HO2.  Deactivation  by  V-T  transfer  was  found  to  be  more  important  for  lower  collision 
temperatures.  The  authors  concluded  vibration-vibration  (V-V)  transfer  was  a  minor 
effect,  because  they  found  negligible  post-collision  excitation  of  the  O2  oscillator. 

From  the  studies  discussed  above  involving  larger  molecules,  there  is  some  evidence 
that  the  characteristics  of  vibrational  relaxation  by  collisions  with  a  bath  gas  change  at 
higher  pressures.  Based  on  the  experimental  results  of  Schwarzer  et  al.  and 
Benzler  et  al,  there  appears  to  be  a  slowing  of  relaxation  time  constants  with  increasing 
bath  density  (pressure).  Similar  density  effects  on  relaxation  time  constants  were  seen  in 
the  trajectory  study  by  Paul,  et  al  We  wish  to  determine  if  these  density  dependent 
effects  are  a  general  phenomenon.  Combustion  conditions  can  occur  at  elevated  densities 
(pressures)  and  it  is  of  interest  to  assess  whether  the  IBC  prediction  of  linear  scaling  of 
relaxation  time  constants  as  a  function  of  system  density  holds  true  in  the  relaxation  of 
small  radicals  and  molecules  important  to  hydrocarbon  combustion.  As  an  initial  system. 
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we  chose  to  study  the  relaxation  of  vibrationally  excited  HO2  in  a  dense  gas  at  a  single 
temperature  with  varying  densities  (pressures).  We  suspect  that  the  single-collision 
relaxation  studies  of  HO2  done  so  far  cannot  offer  conclusive  guidance  about  the 
timescales  of  the  relaxation  process,  or  whether  the  density  of  the  bath  would  have  any 
effect  on  the  relaxation  mechanism.  The  methods  are  outlined  in  Sec.  B,  the  results  in 
Sec.  C,  and  the  conclusions  in  Sec.  D.  Future  studies  are  proposed  in  Sec.  E. 

B.  Methods 

We  have  investigated  the  relaxation  of  vibrationally  excited  HO2  in  a  dense  gas  with 
classical  trajectories.  The  methods  used  in  the  simulations  are  similar  to  the  general 
strategies  developed  by  others  which  have  been  used  to  investigate  the  effect  of  a  dense 
gas  environment  on  isomerization, vibrational  relaxation,'^*^‘^^’'^^  and  unimolecular 
decay. 

38 

All  of  the  simulations  were  performed  using  the  classical  trajectory  code  GENDYN 
modified  for  the  dense  gas  simulations.  The  simulations  of  HO2,  embedded  in  a  bath 
across  a  series  of  pressures  (35,  70,  100,  500,  900,  and  1300  atm),  consisted  of  125  Ar 
atoms  with  a  single  HO2  radical.  The  bath  temperature  was  800  K  and  each  ensemble 
consisted  of  1,600  trajectories.  Additional  tests  with  larger  bath  sizes  of  250,  500,  750 
and  1000  Ar  atoms  were  also  completed  at  1300  atm.  The  velocity  Verlet  integrator  was 
used  for  integrating  the  equations  of  motion.  The  integration  step  size  was  0.1  fs  and  the 
maximum  integration  time  needed  to  relax  HO2  was  determined  from  preliminary 
simulation  results.  The  use  of  the  velocity  Verlet  integrator  was  motivated  by  the  low 
execution  cost  and  long,  stable  integration  behavior  shown  and  discussed  in  Ch.  II  Sec. 
D.  Phase  space  points  of  the  system  were  recorded  every  3  fs  for  analysis. 
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1.  Intermolecular  Potentials 

The  potential  energy  for  the  system  is  given  by 

V  =  V  -\-V  -\-V 

^  ^  HO2  ^  ^  Inter  ^  *  Cell  ’ 

where  is  represented  by  the  XXZLG  and  is 


(23) 


The  heteroatomie  parameters  oiAjj,  By  and  Q  were  obtained  using 


1  T  A 

The  values  for  the  Buekingham  parameters  are  given  in  Table  VII.  Mieheal  et  al. 
eommented  that  the  inclusion  of  electrostatic  forces  would  lead  to  improved  accuracy  in 
modeling  the  intermolecular  interactions  of  HO2  with  the  bath  gas  atoms.  However,  such 
considerations  are  beyond  the  scope  of  the  present  study. 
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The  term  refers  to  the  spherieal  cell  potential  that  was  used  to  confine  the 

131 

simulation.  This  potential  was  developed  by  Marks  et  al.  and  is  a  modified  Bom- 
Mayer  function 

2 

V,ell  =  Z,  ^  ))  (30) 

where  a  is  the  potential  energy  at  the  radius  R  of  the  cell,  n  is  the  distance  of  the  P  atom 
from  the  center  of  the  cell,  and  is  a  parameter  that  controls  the  steepness  of  the  Ar-  and 
HOi-wall  potentials.  The  radius  of  the  cell  varied  with  the  desired  pressure  of  simulation, 
while  «and  y^were  set  to  75  kcal/mol  and  6.5  A' ^  respectively. 

Table  VII.  Values  of  the  intermolecular  potential  parameters. 


0" 

Ar 

A  (kcal/mol) 

2202.082 

69416.3 

69600.0 

B  (A-') 

3.74 

3.96 

4.04 

C  (kcal  mof^  A'^) 

32.5956 

347.3498 

1295.0 

The  parameters  for  H  and  O  were  taken  from  D.  C.  Sorescu,  B.  M.  Rice,  and  D.  L. 
Thompson,  J.  Phys.  Chem.  B  101,  798  (1997). 

The  computational  cost  of  evaluating  the  (n/2)(n-l)+3n  interactions  scales  as  n  , 
where  the  first  term  is  the  number  of  interactions  between  bath  gas  atoms,  the  second 
term  is  the  number  of  interactions  between  the  reactant  atoms  and  bath  atoms,  and  n  is 
the  number  of  bath  gas  atoms.  The  computational  cost  of  the  potential  evaluations  was 
minimized  by  only  evaluating  the  potential  energies  between  atoms  that  had  an 
appreciable  interaction.  This  was  accomplished  by  smoothly  turning  off  the  potential  and 
ignoring  any  interactions  greater  than  a  distance  of  2.7cry  using  a  switching  function  that 
has  two  continuous  derivatives.  The  form  of  this  switching  function  is 
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5(A)  =  min[l, max[ 0,1  -lOgir^f  + 1  Sgir^Y  - 6g{r,.y ]]  (31) 

where 

g{rij)  =  inj  -2.5<Jij){2.1  <Jij  -2.5  (32) 

where  the  Cy  are  the  van  der  Waals  radii  and  g(rij)  is  the  progress  variable  of  the 

switehing  funetion.  The  switehing  funetion  was  only  applied  when  g(rij)  was  between  0 

and  1.  The  value  of  the  unmodified  potential  at  the  start  of  the  attenuation  (r,y=  2.5  Cy)  is 

~3%  of  the  magnitude  of  the  potential  well  depth.  At  the  end  of  the  attenuation  {Xij  = 

2.1  (Jij)  the  unmodified  potential  is  ~2%  of  the  magnitude  of  the  potential  well  depth. 

The  van  der  Waals  radii  of  the  Exp-6  potentials  were  loeated  using  biseetion.  The  use 

of  the  switehing  funetion  should  have  only  minor  effeets  on  the  physieal  properties  of  the 

fiuid.^° 

1  T9 

A  neighbor  list  was  also  used  to  further  reduee  the  eost  of  the  foree  ealeulations  by 
exeluding  pairwise  interaetions  outside  the  eutoff  of  2.1  (jy.  The  organization  of  the 
neighbor  list  is  shown  in  Table  VIII.  The  neighbor  list  eonsists  of  two  arrays:  The  first 
eontains  the  number  of  neighbors  of  a  particular  atom  and  the  second  contains  the 
identification  number  of  the  neighboring  atom. 
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Table  VIII.  A  schematic  of  a  neighbor  list  using  a  two-an  ay  system.® 


®  Based  on  Fig.  5.5  in  Haile,  Ref  30c. 


Updatmg  of  the  neighbor  list  was  automated  monitoring  a  buffer  region  out  to  2.9(7^  and 
defining  a  maximimi  displacement  r„ax  =  2.9<Jij-2.1  <Jy.  When  the  sum  of  the  two  largest 
atomic  displacements,  measmed  fiom  the  most  recent  list  update,  became  greater  than  or 
equal  to  r„ax  llie  list  was  updated. 

Evaluation  of  initial  conditions  and  a  trajectory  of  1  ps  diuatiou,  calculated  for  125  Ar 
atoms  at  800  K  and  1300  atm  with  one  HO2  and  no  reduction  in  the  8,125  pauwise 
interactions  requhes  »310  s  of  cpu  time.  Applying  the  switching  fimction  and  neighbor 
list  in  a  similar  trajectory  reduced  the  average  number  of  interactions  to  =1180,  and  the 
evaluation  time  of  this  tr  ajectory  was  =80  s  of  cpu  time. 

In  oin  previous  work  on  the  isolated  radical  in  Ch.  Ill,  the  gradients  fiom  the  XXZLG 
were  calculated  by  finite  difference  using  a  central  difference  formula.  Also  in  that  work 
we  identified  a  cusp  in  the  potential  along  a  perpendicular  bisector  with  the  origin  bemg 
the  center  of  the  0-0  bond.  In  the  region  of  the  cusp,  we  symmetrized  the  potential  and 
foimd  it  sufficiently  smoothed  the  discontinuity.  Unfortimately,  this  increased  the 
computational  cost  of  the  siuface.  In  an  effort  to  uuprove  the  efficiency  of  oin  simirlatiou 
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in  the  present  work,  the  gradients  of  the  XXZLG  potential  were  ealeulated  analytieally  in 
all  regions  away  from  the  perpendicular  bisector.  In  the  region  of  the  perpendicular 
bisector  we  evaluate  the  gradients  by  finite  difference  on  a  symmetrized  surface  as 
described  in  Ch.  III.  It  was  necessary  to  use  the  mixed  numeric  and  analytic  gradients 
because  of  the  previously  discussed  deficiencies  of  the  XXZLG  PES,  see  Appendix  II. 
Using  analytic  gradients  offers  an  order  of  magnitude  speedup  over  numerical  gradients. 
The  evaluation  times  for  10  ps  trajectories  including  initial  conditions  that  typically 
involved  -675,000  potential  and  gradient  calls  are  as  follows:  356  s  for  analytic  gradients 
with  symmetrized,  numerical  gradients  evaluated  in  the  bisector  region;  233  s  for  pure 
analytic  gradients  without  symmetrization  in  the  bisector  region,  a  procedure  that  yields 
non-energy  conserving  trajectories;  and  2847  s  for  numerical  evaluation  of  the  gradients 
everywhere  including  the  symmetrized  bisector  region.  The  relative  speed  up  from  using 
the  mixed  analytic  and  numeric  gradients  is  a  factor  of  eight,  while  the  full  analytic 
gradients  offer  a  speed  up  factor  of  12.  The  loss  of  efficiency  from  the  mixed  gradient 
solution  is  due  to  the  overhead  of  multiple  potential  energy  evaluations  needed  for 
calculating  both  the  finite-difference  gradients  and  the  symmetrization  of  the  potential. 

2.  Initial  Conditions 

Simulations  of  vibrationally  excited  HO2  in  a  bath  gas  were  performed  for  a  bath 
temperature  of  800  K  and  pressures  of  35,  70,  100,  500,  900  and  1300  atm.  The  initial 
conditions  were  done  in  three  sequential  steps:  (a)  Populate  the  simulation  vessel  with  the 
desired  number  of  bath  gas  atoms,  followed  by  adjustment  of  the  simulation  sphere 
volume  to  obtain  the  desired  pressure  (done  once  for  each  pressure);  (b)  select  the  initial 
conditions  of  the  HO2;  and  (c)  thermalize  the  bath  gas.  After  the  third  step  an  NVE 
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trajectory  was  integrated  using  the  velocity  Verlet  integrator  for  a  predetermined  time 
with  a  step  size  of  0.1  fs.  The  termination  time  of  a  trajectory  was  based  on  approximate 
thermalization  times  from  early  test  simulations.  The  lower  pressure  ensembles  were 
integrated  for  longer  times  to  ensure  sufficient  relaxation  of  the  HO2. 

a.  Cell  Size  and  Pressure 

The  following  paragraphs  are  devoted  to  discussion  of  populating  the  reaction 
environment  and  adjusting  the  volume  to  obtain  a  desired  pressure.  The  initial  volume  of 
the  simulation  vessel  was  calculated  using  the  ideal  gas  law  for  the  desired  temperature 
and  pressure.  This  volume  was  used  as  a  starting  point  for  populating  the  locations  of  the 
Ar  atoms  within  the  confining  potential.  The  simulation  cell  was  centered  at  the  origin  of 
a  Cartesian  coordinate  frame  and  the  center  of  mass  of  the  radical  was  also  positioned 

1  T-5 

there.  The  equilibrium  geometry  of  the  radical  was  obtained  using  the  Stepit  routine 
interfaced  to  GENDYN.  The  gradients  and  Hessians  of  the  PES  were  calculated  at  this 
geometry  to  confirm  this  was  the  minimized  structure  of  the  radical. 

The  bath  gas  atoms  were  randomly  inserted  in  the  simulation  cell  using  a  simple 
acceptance/rejection  scheme.  A  trial  placement  of  a  bath  gas  atom  was  done  by 
randomly  selecting  the  Cartesian  coordinates  of  one  of  the  previously  placed  atoms, 
adding  a  random  displacement  to  all  three  Cartesian  coordinates,  inserting  a  test  atom  at 
the  resulting  position  and  evaluating  of  the  potential  energy.  The  maximum  magnitude  of 
random  displacement  in  a  single  direction  was  taken  to  be  half  of  the  current  cell  radius. 
The  new  placement  was  accepted  if  the  total  potential  energy  remained  less  than  zero; 
otherwise  the  placement  was  rejected  and  a  new  trial  placement  was  generated.  This 
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process  was  repeated  until  all  of  the  bath  gas  atoms  were  placed  in  the  simulation  cell.  A 
larger  number  of  insertion  rejeetions  oecurred  for  the  higher  pressure  simulations,  whieh 
required  smaller  simulation  volumes.  A  smaller  vessel  volume  has  two  effects  on  the 
insertion  proeedure:  First,  there  was  an  increased  ehanee  of  sampling  of  the  repulsive 
spherieal  wall;  second,  the  sequential  plaeement  of  bath  gas  atoms  eventually  leads  to 
paeking  of  the  bath  gas  in  which  insertions  with  positive  potential  energies  are  more 
likely.  A  simple  solution  to  prevent  large  rejection  rates  was  to  start  with  a  larger 
estimated  ideal  gas  volume  and  eorrect  for  the  overestimated  volume  in  the  pressure 
adjustment  trajeetory.  The  lowest  aeeeptance  rate  for  the  eell  population  proeedure  was 
-20%. 


After  placing  all  of  the  bath  gas  atoms,  a  trajectory  was  integrated  for  25  ps  with  a  0.1 
fs  step  size  and  the  pressure  was  equilibrated  by  changing  the  volume  of  the  simulation 
cell.  The  HO2  was  held  stationary  at  the  eenter  of  the  cell  during  this  proeess.  This 
pressure  equilibration  trajectory  was  started  by  assigning  initial  Cartesian  veloeities  to  the 
bath  gas  from  one-dimensional  Maxwell  veloeity  distributions.  Maxwell  distributions 
were  calculated  by  sampling  a  weighted  Gaussian  distribution.  The  weighting  faetor  used 
is  (miKT)^^^  where  m,  is  the  mass  of  atom,  k  is  Boltzmann’s  constant,  and  T  is  the 
temperature.  The  Gaussian  distributions,  for  eaeh  Cartesian  degree  of  freedom,  were 
calculated  using  the  gasdev  subroutine  from  Numerical  Recipes.  Velocities  were 
assigned  to  the  bath  gas  while  the  radical  was  held  fixed  at  the  eenter  of  the  simulation 
cell  by  exeluding  it  from  the  integration  of  the  equations  of  motion;  however,  the 
pairwise  interactions  of  the  radical  with  the  Ar  atoms  were  ealculated  in  the  equations  of 
motion  to  prevent  any  aphysical  overlaps. 
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After  selecting  the  initial  velocities  of  the  bath  atoms,  a  trajectory  was  integrated  and 
after  each  integration  step  the  instantaneous  average  kinetic  energy  (temperature)  was 
calculated  using^® 


=— E 


2  ,  2  ,  2 
Px  +Py+P. 


(33) 


SNk  —  m 

where  N  is  the  number  of  bath  atoms,  and  px,  Py,  Pz  are  the  Cartesian  momentum 
components  of  each  atom.  After  evaluating  the  temperature,  the  pressure  was  calculated 

•  30 

using 


where  p  is  the  current  number  density,  7)  is  the  result  of  Eq.  33,  E  is  the  current  volume 
of  the  sphere  at  the  value  of  R  in  equation  30,  and  Y!E/ij'P(^ij)  is  the  virial  from  the 
interatomic  forces  F(rij).  Both  the  temperature  and  pressure  were  smoothed  by  a  rolling 
average^^"^  with  a  200  fs  window  (2,000  integration  steps): 

=  -E" +  G(n  - 1)  + ...  +  G{i)  -  G{i  - 1) ,  (35) 

n 

where  each  G(i)  is  either  an  instantaneous  temperature  or  pressure,  n  is  the  total  number 
of  samples  in  the  window  and  Groii  is  the  value  of  the  rolling  average  of  either 
temperature  or  pressure.  After  every  1000  integration  steps  (100  fs),  the  rolling  average 
temperature  of  the  bath  was  compared  to  the  desired  temperature  and  if  the  absolute 
difference  was  greater  than  10  K  the  momenta  were  scaled  using 
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where  and  Pq  are  the  new  and  old  momenta,  respectively,  and  T  is  the  desired 
temperature.  After  the  temperature  was  tested  the  rolling-average  pressure  was  compared 
to  the  desired  pressure.  If  the  absolute  difference  was  greater  than  5%  of  the  desired 
pressure  then  the  volume  of  the  reaction  cell  was  resized  to  give  the  desired  pressure. 
The  Cartesian  positions  of  the  bath  gas  atoms  were  scaled  to  adjust  for  the  change  in 
volume  using 


(In  =(1o 


ry 

yVoj 


(37) 


where  qN  and  qo  are  the  new  and  old  Cartesian  coordinates,  respectively,  and  and  Vq 
are  the  new  and  old  cell  volumes.  The  cell  radius  was  smoothed  with  a  50  point  rolling 
average  using  Eq.  37,  as  a  function  of  the  number  of  times  it  was  adjusted.  As  stated 
above  the  cell  volume  was  only  adjusted  if  the  absolute  difference  between  the  rolling 
average  and  desired  pressure  differed  by  5%. 


Figure  14  shows  examples  of  the  instantaneous  (magenta  curves)  and  rolling  average 
values  (black  curves)  for  cell  radius  (a),  temperature  (b),  and  pressure  (c)  as  functions  of 
time.  The  temperature  and  pressure  targets  are  800  K  and  1,300  atm.  The  instantaneous 
and  rolling  average  values  oscillate  about  the  target  values  with  the  magnitude  of  the 
rolling  average  values  being  modestly  smaller  than  the  instantaneous  values  for  both 
temperature  and  pressure.  The  rolling  average  of  the  cell  radius  reaches  a  smooth 
average  within  1  ps.  Frames  (d-f)  are  an  enhanced  view  on  a  shorter  timescale  of  the 
results  shown  in  (a-c).  These  plots  (d-f)  show  that  the  large  initial  fluctuations  in  the 
temperature,  pressure,  and  the  cell  radius  are  quickly  smoothed  and  stable  trends  for  all 
three  properties  are  reached  well  within  the  maximum  integration  time  of  25  ps.  The  cell 
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radius  was  fixed  to  the  final  value  obtained  at  the  end  of  the  integration  time.  Similar 
results  for  the  smoothing  of  the  cell  fluctuations,  as  those  shown  in  Fig.  14,  were  obtained 
for  all  of  the  temperatures  and  pressures  simulated  in  this  work.  The  averaged  volumes, 
radii,  and  bath  densities  are  reported  in  Table  IX. 


(d)  - 

A 

-  V  ^ 

■  ■ 

Figure  14  An  illustration  of  the  convergence  of  the  cell  properties  during  the  pressure 
equilibration;  (a)  Cell  radius;  (b)  Bath  temperature;  (c)  Bath  pressure.  On  a  shorter 
timescale:  (d)  Cell  radius  (e)  Bath  temperature;  (f)  Bath  pressure 

Table  IX,  The  radii,  volumes,  and  gas  densities  of  the  simulated  pressures.^ 


Pressure  (atm) 

Radius  (A) 

Volume  (A^) 

Density  (mol/L) 

35 

45.62  (4)" 

398000  (1200) 

0.522  (2) 

70 

36.26  (7) 

120000  (1200) 

1.039(6) 

100 

32.34  (6) 

141800 (700) 

1.46(1) 

500 

19.67  (9) 

31900  (400) 

6.51  (9) 

900 

16.78  (7) 

19800  (300) 

10.5  (1) 

1300 

15.36(7) 

15200  (200) 

13.7(2) 

These  values  are  the  averages  and  uncertainties  calculated  from  the  sixteen  sub¬ 
ensembles. 

The  numbers  in  parenthesis  indicate  the  uncertainty  in  the  last  non-zero  digit. 

The  volume  of  the  cell  decreased  by  a  factor  of  three  going  from  35  to  1300  atm  and  the 

density  increased  two  orders  of  magnitude.  At  the  higher  pressures  of  900  and  1300  atm 
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there  appears  to  be  some  signs  of  incompressibility  with  only  modest  changes  in  cell 
radius  and  density. 

b.  Solute  Initial  Conditions 

The  initial  conditions  of  the  radical  were  selected  using  with  a  total 

internal  energy  of  40  kcal/mol  and  no  restriction  on  the  angular  momentum  of  the  radical. 
The  bath  atoms  were  held  fixed  at  the  positions  reached  at  the  end  of  step  (a)  and  only  the 
reactant  was  changed  during  the  Markov  walk.  The  intermolecular  potential  was  ignored, 
so  the  surrounding  bath  gas  had  no  effect  on  the  selection  of  the  internal  state  of  the 
reactant.  Any  possible  interatomic  overlaps  between  HO2  and  Ar  that  could  have 
occurred  were  dealt  with  in  step  (c).  A  warm-up  walk  of  250,000  steps  in  Cartesian 
coordinates  was  started  from  the  equilibrium  geometry  in  which  all  three  atoms  of  HO2 
were  allowed  to  move  with  random  displacements.  The  maximum  possible  displacement 
of  a  single  atom  was  0.2  A.  The  acceptance  of  a  new  displacement  was  governed  by  the 
weight  in  Eq.  7.  At  the  end  of  the  walk  a  set  of  momenta  were  sampled  and  scaled  to 
obtain  the  desired  energy  as  described  in  Ch.  II.  The  final  coordinates  of  HO2  at  the  end 
of  the  Markov  walk  were  used  as  the  starting  configuration  for  selecting  the  initial 
conditions  of  the  next  trajectory  in  the  simulation.  Subsequent  walk  lengths  consisted  of 
100,000  steps  that  were  taken  between  trajectories.  The  length  of  subsequent  walks  was 
selected  based  on  the  convergence  of  the  coordinates  and  potential  energy  as  a  function 
of  the  number  of  steps  in  the  Markov  walk.  The  number  of  accepted  moves  of  each  walk 
was  between  0.4  and  0.6. 


102 


119 


Assigning  the  translational  momenta  is  the  next  step  in  selecting  the  initial  conditions 
of  the  solute.  A  set  of  Maxwell  velocity  distributions  were  sampled  for  each  of  the  atoms 
in  HO2,  then  the  center-of-mass  velocity  was  calculated  using 

(38) 

where  pi  are  the  Cartesian  momenta  and  mi  is  the  mass  of  the  ith  atom.  The  components 
of  the  sampled  velocities  that  contribute  to  the  internal  motions  of  HO2  were  discarded. 
Translational  momentum  was  then  vectorially  added  to  the  momenta  sampled  at  the  end 
of  the  EMS  procedure  using 

Pi=Pi-mv^.  (39) 

The  result  is  a  rotationally  and  vibrationally  hot  reactant  that  has  a  thermal  distribution  of 

translational  energy.  The  sampled  kinetic  energy  of  the  center-of-mass  motion  was 
compared  to  the  analytic  Maxwell-Boltzmann  distribution  of  kinetic  energy  given  by 

P(EJ  =  2j - A^e?qi  - ^  ,  (40) 

\7:(KTf  yKT) 

where  Et  is  the  center-of-mass  kinetic  energy.  An  example  histogram  of  the  translation 
kinetic  energy  distribution  is  shown  in  Fig.  15. 
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Figure  15  The  histogrammed  center-of-mass  kinetic  energy  of  HO2  from  the  selection  of 
the  initial  conditions  at  800  K  shown  as  red  boxes.  The  analytic  distribution  for  T  =  800 
K  is  shown  as  a  solid  black  line  for  comparison. 

c.  Final  Thermalization 

The  conformation  of  the  radical  was  changed  by  the  selection  of  microcanonical 
initial  conditions  described  in  Sec.  B.2b,  without  consideration  of  the  configuration  of  the 
surrounding  bath  gas.  Interatomic  overlaps  between  the  HO2  and  Ar  atoms  are  a  possible 
consequence  that  could  occur  due  to  this  neglect.  To  correct  for  this  possibility,  the  bath 
gas  atoms  were  re-thermalized  with  a  new  set  of  velocities  using  the  procedure  as 
described  in  Sec.  B.2a.  The  reactant  was  again  held  fixed,  as  described  in  Sec.  B.2a  but 
retaining  the  phase  space  point  selected  during  Sec.  B.2b.  A  thermal  trajectory  was 
integrated  for  25  ps  with  the  temperature  controlled  in  the  same  manner  as  it  was  in  Sec. 
B.2a.  The  final  Cartesian  coordinates  of  the  bath-gas  atoms  were  stored  at  the  end  of 
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each  thermalization  trajectory,  and  used  as  the  starting  configuration  for  subsequent 
thermalization  trajectories  in  the  simulation  sequence. 

3,  Energy  Transfer 

Analysis  of  the  energy  transfer  between  the  hydrogen  peroxyl  radical  and  the  bath 
was  done  by  decomposing  the  Hamiltonian  of  the  radical  into  the  translational,  rotational, 
and  vibrational  components.  The  separation  of  linear  momenta  of  the  radical  is  exact. 
After  removing  the  translational  component  a  zeroth-order  separation  of  rotational  and 
vibrational  energy  was  assumed 

H,,=H,+H,,  (41) 

where  Hr  the  was  calculated  using 

(42) 

L  is  the  whole-body  rotational  angular  momentum 

L  =  rx  p,  (43) 

and  the  instantaneous  moment  of  inertia  was  calculated  using 

m{y^+z^)  -mxy  -mxz 

/=  -mxy  m{x^+z^)  -myz  .  (44) 

-  mxz  -  myz  m(x^  +  ) 

This  allows  the  vibrational  energy  to  be  defined  as 
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(45) 

Frederick  et  al.  described  several  different  ways  of  separating  the  vibrational  and 
rotational  degrees  of  freedom  in  molecular  systems  and  recommended  using  Eq.  42, 
because  the  definition  is  independent  of  the  coordinate  system. 

The  vibrational  energy  was  calculated  using  Eqns.  42-45  as  a  function  of  time  and  the 
asymptotic  thermal  energy  E'(oo)  was  subtracted.  We  refer  to  the  vibrational  energy 
above  E'(oo)  simply  as  “excess-vibrational  energy.”  We  assume  that  two  factors  will 
govern  the  behavior  of  E'(oo)  in  the  limit  of  infinite  time:  Eirst,  the  solute  will  have  the 
same  average  thermal  energy  as  the  bath.  Second,  in  this  limit  the  radical  will  have  lost 
sufficient  energy  such  that  the  trajectory  is  only  sampling  harmonic  portions  of  the  PES 
and  that  the  contribution  of  the  potential  to  equipartion^^^  will  be  ^/dT.  The  definition  of 

E{<x>)  is  3kT  for  vibrational  motion  and  fx"!  for  rotational  motion.  Since  we  are 
simulating  a  finite  system,  the  energy  loss  from  the  radical  will  cause  the  bath  to  heat  to  a 
new  asymptotic  temperature;  consequently,  we  assigned  T  to  be  the  time-dependent 
ensemble  average  temperature  of  the  bath  gas. 

The  resultant  curve,  representing  the  decay  of  the  excess-vibrational  energy  of  the 
radical  to  E{<x)),  was  fit  to  a  bi-exponential  function 

E{t)  =  A  exp  +  B  exp  {—kA) ,  (46) 

where  A  and  B  are  calculated  by  linear  least  squares  subject  to  the  constraint  A+B  =  1, 

T9 

and  ki  and  k2  are  time  constants  fitted  using  standard  nonlinear  methods. 
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C.  Results 

1.  Thermal  State  of  Bath 

The  first  test  of  the  simulation  method  was  to  show  that  the  equilibration  of  the  bath 
resulted  in  a  stable  temperature  and  pressure  and  that  the  Maxwell-Boltzman  distribution 
of  speeds  was  adequately  sampled  by  the  bath  gas.  These  tests  were  done  to  ensure  that 
there  were  no  abnormalities  in  the  implementation,  sinee  previously  Marks  et  al.  had 
shown  that  satisfactory  results  for  these  properties  can  be  obtained  using  the  sphere 
potential  with  a  smaller  number  of  bath  atoms.  We  tested  two  temperatures  above  and 
below  800  K  for  additional  verification.  The  initial  conditions  of  HO2  were  selected  so 
that  10  kcal/mol  was  divided  between  the  rotational  and  the  vibrational  degrees  of 
freedom  and  the  translational  velocities  were  selected  from  a  thermal  distribution,  as 
described  in  Sec.  B.2.  The  magnitude  of  the  internal  energy  was  of  HO2  was  arbitrarily 
selected  to  be  smaller  than  the  value  (40  kcal/mol)  that  used  in  the  actual  relaxation 
studies.  This  choice  was  made  to  minimize  potential  heating  of  the  bath  gas  from  energy 
transfer  out  of  HO2.  It  should  also  be  pointed  out  that  the  results  reported  in  this  section 
are  from  individual  trajectories. 

Figure  16  (a-c)  contains  of  histograms  of  the  Maxwell-Boltzmann  speed  distribution 
calculated  from  the  bath  gas  speeds  sampled  during  a  trajectory,  the  analytic  function  is 
also  plotted  as  dashed  lines  for  comparison.  These  plots  were  generated  from  the  results 
of  individual  trajectories  of  125  Ar  atoms  and  a  single  HO2. 
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Figure  16  Thermal  state  and  pressure  of  the  bath  for  a  variety  of  sampled  temperatures, 
(a)  Maxwell-Boltzmann  bath  gas  speed  distribution;  (a)  400  K;  (b)  800  K;  (e)  1600  K. 
Instantaneous  pressures;  (d)  35  atm;  (e)  1300  atm;  (f)  500  atm.  Instantaneous 
temperatures;  (g)  400  K,  (h)  800  K,  (i)  1600  K. 

The  velocity  histograms  are  in  excellent  agreement  with  the  analytic  distributions. 
The  temperature  and  pressure  calculated  using  Eqns.  33  and  34  are  shown  in  frames  (g-i) 
and  (d-f),  respectively.  The  fluctuations  in  the  pressure  are  a  function  of  both  the 
temperature  and  the  magnitude  of  the  interatomic  forces  (see  Eq.  34),  so  as  the 
temperature  and  density  of  the  system  increase  the  magnitude  of  the  pressure  fluctuations 
also  increases.  However,  it  is  either  the  time  average  or  ensemble  average  that  has 
physical  meaning.  The  time  average  pressures  going  down  the  middle  column  of  Eig.  16 
are;  34.8  atm,  1242.1  atm,  and  506.6  atm  and  the  root  mean  square  (rms)  fluctuations  are; 
1.8,  130.4,  and  26.3  atm,  respectively.  The  time  average  temperatures  in  the  right  hand 
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column  of  Fig.  16  are:  396.8  K,  111 .1  K,  and  1615.8  K  and  the  rms  fluctuations  are  6.0 
K,  20.0  K,  and  23.4  K,  respeetively.  The  increased  fluctuations,  in  the  temperatures  and 
pressures,  at  higher  densities  eould  be  indieative  of  increased  incompressibility  of  the 
bath  gas.  This  incompressibility  can  be  identified  as  an  inereased  number  of  multiple, 
non-negligible  interatomie  interactions.  Some  evidence  of  this  incompressibility  can  be 
seen  in  the  rms  fluctutations  of  the  HOi-Ar  interactions  calculated  by  Eq.  25:  The  rms 
fluctuations  in  the  HOi-Ar  interaetions  are  the  following:  for  400  K  and  35  atm,  it  is  0.16 
keahmol;  for  800  K  and  1300  atm,  it  is  0.90  keal/mol;  and  for  1600  K  and  500  atm,  the 
rms  is  0.79  kcal/mol.  The  multiple,  non-negligible  interactions  cause  more  fluctuation  in 
the  kinetie  energy  of  the  system,  which  is  used  to  define  a  temperature  in  Eq.  33.  The 
virial  term,  used  in  caleulating  the  system  pressure  in  Eq.  34,  also  fluctuates 

more,  which  combined  with  the  modulation  in  temperature  could  be  the  eause  the 
observed  fluctuation  in  pressure. 

There  is  a  slight  upward  drift  in  the  temperature  of  the  bath  for  the  800  K  (frame  h) 
and  1600  K  (frame  i).  The  drift  is  due  to  the  bath  gas  atoms  gaining  energy  from  HO2  via 
collisions  during  the  course  of  the  trajectory.  The  expected  thermal  energy  of  rotational 
and  vibrational  energy  of  HO2  is  4.5A:r.  Eor  the  800  K  trajeetory  the  temperature  of 
drops  as  low  as  =720  K  and  the  relative  internal  energy,  10  kcahmol,  of  HO2  at  this 
temperature  is  ~6.9kT,  so  it  is  probable  that  energy  will  transfer  to  the  bath  gas  atoms  and 
cause  the  overall  temperature  of  the  simulation  to  rise.  Similar  analysis  1600  K  trajectory 
shows  the  same  behavior,  the  instantaneous  temperature  drops  as  low  as  =1510  K  and  the 
relative  internal  energy  of  HO2  is  ~4.6kT,  so  again  energy  transfer  out  of  HO2  ean  oceur. 
The  400  K  trajectory  shows  no  overall  temperature  drift  beeause  the  aetivating  and 
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deactivating  collisions  of  HO2  with  the  bath  gas  atoms  effectively  canceled  each  other  out 
over  the  course  of  the  trajectory. 

A  way  to  decrease  the  magnitude  of  the  fluctuations  is  to  increase  the  number  of  bath 
gas  atoms.  If  one  assumes  that  the  reduction  in  error  follows  the  relation  of  1/-//V, 
where  N  is  the  system  size,  reducing  the  fluctuations  in  Fig  16  by  half  would  require  a 
system  size  of  500  Ar  and  reducing  the  error  to  a  quarter  of  its  value  would  require  2000 
Ar.  These  results  taken  as  a  whole  show  that  the  procedure  described  in  Sec.  B.2b  was 
sufficient  in  setting  the  thermodynamic  state  of  the  bath. 

2,  Energy  Transfer  Results 

The  goal  of  this  study  is  to  detemine  whether  there  is  any  pressure  effect  on  the 
relaxation  of  excited  HO2.  The  integration  time  lengths  were  determined  from 
extrapolation  of  preliminary  results.  The  35  and  70  atm  trajectory  ensembles  were 
simulated  for  1700  ps,  the  100  atm  ensemble  was  simulated  for  750  ps,  and  the  500,  900, 
and  1300  atm  ensembles  were  simulated  for  500  ps.  The  effect  of  the  system  size  was 
studied  by  varying  numbers  of  bath  gas  atoms  in  the  simulation  vessel.  The  number  of 
bath  gas  atoms  was  varied  from  125  up  to  1000  Ar  atoms  and  these  simulations  were 
done  at  1300  atm  and  800  K. 

a.  Vibrational  Relaxation  Pressure  Dependence 

Figure  17  contains  plots  of  the  vibrational  energy  of  HO2,  the  bath  temperature,  and 
the  bath  pressure  as  functions  of  time.  All  of  these  quantities  were  scaled  to  a  range  of  0 
to  1  relative  to  their  respective  values  at  time  =  0.  The  initial  conditions  of  HO2  were 
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selected  so  that  the  translational  motion  was  thermal  and  40  kcal/mol  was  distributed 
between  the  rotational  and  vibrational  degrees  of  freedom.  Both  the  bath  temperature  and 
pressure,  for  all  of  the  trajectory  ensembles,  show  a  relative  increase  with  respective  to 
the  values  at  time  =  0.  The  final  relative  increase  in  both  temperature  and  pressure,  for 
all  of  the  ensembles,  is  ~10%.  The  increase  in  both  temperature  and  pressure  is  due  to 
the  Ar  atoms  gaining  energy  as  the  HO2  is  deactivated  by  collisions. 
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Figure  17  ,  The  results  of  HO2  with  40  kcal/mol  of  excess  internal  energy  relaxing  in  a 
bath  of  125  Ar  atoms  at  various  pressures.  The  scaled;  (a)  vibrational  energy,  (b) 
temperature,  and  (c)  pressure.  The  line  coloration;  35  atm,  red;  70  atm,  magenta;  100 
atm,  green;  500  atm,  blue;  900  atm,  yellow;  1300  atm,  dark-green. 
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Figure  18  contains  plots  of  the  energy  decay  curves  and  single  exponential  fits  on  a 
semi-log  scale.  The  amplitude  of  the  single  exponential  fits  was  constrained  to  unity.  It  is 
apparent  that  the  vibrational  relaxation  is  not  well  fit  by  a  single  exponential  and  two 
relaxation  timescales  are  present.  Close  examination  of  the  curves  show  indicate  that 
there  are  two  appearant  linear  regions.  An  example  of  the  transition  from  the  first 
relaxation  timescale  to  the  second  can  be  seen  with  the  red  curve  at  ~  425  ps.  This 
nonlinearity  in  the  energy  relaxation  curves  exists  in  all  of  the  simulated  results. 


Time  (ps) 


Figure  18.  Plots  of  the  decay  of  excess-vibrational  energy  above  the  thermal  limit  (See, 
Fig.  17a)  as  a  function  of  time  on  a  semi-log  scale:  single  exponential  fits  (dashed-black); 
35  atm,  red;  70  atm,  purple;  100  atm,  light-green;  500  atm,  blue;  900  atm,  yellow;  and 
1300  atm  dark-green. 
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The  results  of  fitting  to  Eq.  46  are  shown  in  Fig.  19  and  it  appears  that  the  bi¬ 
exponential  fit  is  satisfaetory.  HO2  loses  at  least  70%  of  the  initial  vibrational  energy  for 
all  of  the  pressures  simulated  and  for  the  highest  pressure  the  vibrational  energy  falls  by 
~95%  within  the  simulated  time. 


Figure  19.  Same  as  figure  18,  exeept  seale  is  not  semi-log  and  fitted  eurves  are  from  Eq. 
46. 

The  fitting  parameters  and  uncertainties  calculated  in  fitting  the  curves  in  Fig.  17a  to 
Eq.  46  are  tablulated  in  Table  X.  The  asymptotic  limit  of  zero  required  for  fitting  the 
curves  in  Fig.  17a  was  done  by  subtracting  off  the  thermal  energy  of  the  radical,  which  is 
a  function  of  the  bath  temperature  shown  in  frame  17(b),  with  an  appropriate  conversion 
by  equipartion,  described  in  methods  Sec.  B.3.  As  discussed  in  the  methods  (Sec.  B.3), 
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the  fitting  of  the  energy  decay  curves  was  done  with  A  and  B  determined  by  linear  least- 
squares  with  the  constraint  of  A  +  B  =  1  and  ki  and  k2  treated  as  adjustable  parameters. 
The  error  bars  in  Table  X  were  generated  by  bootstrapping  the  1600  trajectories  of  each 
ensemble.  In  total,  5000  resampled  groups  of  400  trajectories  were  generated  for  each 
pressure  ensemble.  The  excess-vibrational  energy  curves,  calculated  from  the  resampled 
data  sets,  were  also  fit  to  Eq.  46.  The  bootstrapping  was  done  so  that  each  trajectory 
could  only  be  sampled  once  per  resampled  group. 


Table  X,  Parameters  and  error  bars  from  fitting  the  vibrational  relaxation  of  HO2  to  Eq. 
46 


P  (atm) 

p  (mol/L) 

A 

B 

Oab 

ki^ 

Okl 

k2 

Ok2 

35 

0.52 

0.34 

0.66 

0.0008 

0.0007 

Mjiiijn 

70 

1.04 

0.34 

0.66 

0.07 

0.0010 

100 

1.46 

0.23 

0.77 

0.07 

0.0017 

0.0013 

500 

6.51 

0.21 

0.79 

0.06 

0.0049 

0.0032 

900 

10.47 

0.26 

0.74 

0.06 

0.0051 

0.0045 

1300 

13.66 

0.26 

0.74 

0.05 

BBhIH 

0.0075 

0.0058 

a.  The  unit  of  k],  k2  and  the  respective  errors  is  ps  .  The  amplitudes  A  and  B  are 
dimensionless. 


Eigure  20a  contains  plots  of  the  amplitudes  and  time  constants  found  in  Table  X  as  a 
function  of  the  bath  gas  density.  The  change  of  the  independent  variable  to  density  was 
made  to  faciliate  comparison  with  the  results  of  Paul  et  Both  time  constants  grow 
with  increasing  bath  density,  but  there  is  a  change  in  behavior  near  ~l-2  mol/L  as  shown 
by  the  blue  curves.  Also  at  this  density  the  amplitudes  of  the  exponentials  shown  as  the 
red  curves  oscillate  and  the  amplitude  of  the  primary  exponential  increases,  while  the 
amplitude  of  the  smaller  exponential  decreases.  However,  the  error  bars  of  the  amplitude 
show  a  fair  amount  of  overlap  across  the  density  range  indicating  that  the  amplitudes 
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have  a  weak  density  dependenee.  The  amplitude  of  the  slower  exponential  is  consistently 
larger  than  that  of  the  faster  exponential. 

The  vibrational  relaxation  results  of  Paul  et  are  included  in  frame  (b)  of  Fig.  20 
for  comparison  with  our  results  in  frame  (a).  Plots  of  this  sort  were  not  reported  by  Paul 
et  al,  so  we  have  plotted  their  data  in  Table  III  of  Ref  19c  as  a  function  of  the  reported 
bath  densities.  The  reported  results  are  from  simulations  of  vibrationally  excited  CeFe 
embedded  in  an  N2  fluid  described  with  periodic  boundary  conditions.  Paul  et  al.  fit  the 
relaxation  of  excess-vibrational  energy  was  fit  to  the  bi-exponential  function  where  E('x^) 
was  treated  as  an  adjustable  fitting  parameter.  These  simulations  also  span  a  density 
range  four  times  the  maximum  density  we  investigated.  The  magnitude  of  the  relaxation 
time  constants  for  CeFe  are  larger  than  those  reported  by  us  for  HO2  by  about  a  factor  of 
three  to  five.  We  believe  this  difference  could  be  due  to  CeFe  having  only  low  frequency 
modes  that  will  readily  couple  to  the  motion  of  the  bath  gas.  Despite  the  difference  in  the 
magnitudes  of  the  time  constants  there  are  striking  similarities  in  the  relaxation  results  of 
HO2  and  CeFe.  The  amplitudes  shown  in  Fig  20b  exhibit  a  density  dependence,  but  in  the 
same  range  of  densities  we  simulated,  the  amplitudes  reported  appear  to  exhibit  only  a 
weak-density  dependence.  It  is  interesting  that  the  relaxation  results  for  very  different 
chemical  species  are  well  fit  by  a  bi-exponential  function.  While  there  are  differences  in 
the  magnitudes  of  the  relaxation  timescales,  both  the  CeFe  and  HO2  results  show  two 
different  regions  of  linear  dependence  with  respect  to  bath  density  The  turnover  from 
one  linear  region  to  the  other  in  our  results  occurs  at  ~1  mol/L  and  at  ~10  mol/L  for  the 
results  of  Paul  et  al.  It  is  possible  that  the  size  of  the  molecule  plays  a  role  in  determining 
the  density  at  which  this  turnover  behavior  occurs.  The  turnover  effect  is  somewhat 
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appearant  with  the  dashed  green  line  in  frame  (b),  but  it  is  less  elear  when  examining  the 
solid  green  eurve.  It  is  worth  noting  that  a  turnover  ean  be  seen  in  the  experimental 
relaxation  time  constants  obtained  by  Schwarzer  et  al.  and  Benzler  et  and  in  both 
of  these  studies  the  turnover  occurred  at  a  density  of  ~10  mol/L.  Additional  simulation 
results  for  CeFe  at  a  lower  density  range  are  needed  to  better  characterize  the  density  at 
which  the  turnover  behavior  occurs. 
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Figure  20,  Plots  of  the  fitting  amplitudes  and  time  constants  in  Table  X  as  a  function  of 
bath  gas  density.  The  lines  are  intended  to  guide  the  eye:  (a)  the  solid-blue  line 
corresponds  to  the  larger  time  constant  and  the  amplitude  is  matched  with  the  solid-red 
line.  The  dashed-blue  line  corresponds  to  the  smaller  time  constant  and  the  amplitude  is 
matched  with  the  dashed-red  line;  (b)  the  results  from  Paul  et  al.  (a)  the  solid-green 
line  corresponds  to  the  larger  time  constant  and  the  amplitude  is  matched  with  the 
solid-yellow  line.  The  dashed-green  line  corresponds  to  the  smaller  time  constant  and  the 
amplitude  is  matched  with  the  dashed-yellow  line.  The  symbols  in  (b)  were  included  to 
emphasize  the  densities  specific  densities  of  the  simulations. 
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Our  low-pressure  time  constants  are  shown  on  a  much  larger  scale  in  Fig.  21.  Linear 
fits  are  shown  as  the  black  lines  in  this  figure.  These  fits  to  the  three  data  points  were 
obtained  with  the  intercept  fixed  at  zero.  The  zero  intercept  is  the  limit  of  no  collisions 
and  consequently  no  energy  transfer.  The  fitted  lines  lie  outside  of  the  bootstrapped  error 
bars  reported  at  the  lowest  density.  An  overlap  of  linear  fits  with  the  data  would  have 
been  supporting  evidence  of  a  linear  proportionality  between  the  time  constants  and  the 
bath  gas  density,  which  is  an  important  assumption  in  IBC.  We  think  that  the  increase  in 
both  temperature  and  pressure  of  the  simulation  drives  the  relaxation  of  HO2  and  leads  to 
an  overestimation  of  the  time  constants.  This  possible  finite-size  effect  motivated  a  study 
of  how  the  system  size  affected  the  relaxation  process. 
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Figure  21  The  low-pressure  time  constants  ki  and  k2  as  a  function  of  pressure  on  a 
smaller  scale  than  that  shown  in  Fig.  19.  The  larger  time  constants  are  plotted  relative  to 
the  axis  on  the  left  and  the  smaller  time  constants  are  plotted  relative  to  the  axis  on  the 
right.  The  solid  and  dashed-black  lines  correspond  to  a  linear  fit  of  the  time  constants 
subject  to  the  constraint  of  zero  energy  transfer  in  the  limit  of  zero  pressure. 

b.  Finite  Size  Effects 

The  simplest  prescription  to  avoid  the  finite  size  problem  is  to  increase  the  system 
size  until  the  results  of  interest  no  longer  change.  However,  increasing  the  bath  comes 
with  a  tradeoff  of  time,  because  the  computational  cost  of  the  simulation  scales  with  the 
system  size.  For  our  particular  study,  we  assume  that  the  excited  HO2,  started  with  40 
kcahmol  of  internal  energy,  is  relaxing  to  thermal  equilibrium.  In  the  long-time  limit,  the 
internal  energy  of  the  radical  predicted  by  equipartition  would  be  4.5A:r  or  7.2  kcal/mol  at 
800  K.  This  means  that  ~33  kcal/mol  of  excess  energy  would  be  transferred  to  the  bath 
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in  the  limit  of  complete  thermalization.  The  predicted  increase  in  the  bath  temperature 
for  a  monoatomic  gas  is  -INkIZEu,  so  doubling  the  number  of  bath  atoms  would  reduce 
the  rise  in  temperature  by  half.  To  study  the  effect  of  the  number  of  Ar  atoms  on  the 
relaxation,  we  did  additional  simulations  of  HO2  relaxing  with  increasing  numbers  of 
bath  gas  atoms  with  initial  conditions  selected  described  in  Sec.  B.2b.  at  800  K  and  1300 
atm.  We  chose  to  simulate  several  system  sizes  to  see  how  the  temperature  increase 
progressed:  for  125  atoms  the  increase  would  be  ~90  K;  for  250  atoms,  ~45K,  for  500 
atoms,  ~22  K,  for  750  atoms,  ~15  K,  and  for  1000  bath  atoms,  =^11  K.  We  attempted  to 
partially  account  for  heating  of  the  bath  on  the  relaxation  process  by  making  the 
asymptotic  energy  of  HO2  a  function  of  the  bath  temperature.  However,  this  does  not 
account  for  the  associated  rise  in  pressure  that  occurs  with  heating  (see.  Fig.  17c). 

The  results  from  the  simulations  with  increasing  numbers  of  Ar  atoms  are  shown  in 
Fig.  22.  The  approximate  temperature  increase  shown  in  frame  (a)  agrees  fairly  well  with 
the  predicted  increase  for  each  of  the  numbers  of  bath  gas  atoms.  The  increase  in  the 
simulation  temperature  from  smallest  to  largest  number  of  bath  gas  atoms  is:  36,  18,  12 
and  9  K.  The  increases  in  pressure  shown  in  (b)  are  a  possible  indication  as  to  why  the 
smaller  simulations  have  larger  relaxation  time  constants.  The  pressure  calculated  by  Eq. 
34  has  explicit  temperature  dependence  in  the  ideal  gas  term  pKl.  The  increase  in 
temperature  would  cause  a  corresponding  increase  in  pressure,  or  the  collision  frequency 
between  the  radical  and  the  bath,  which  would  likely  manifest  as  a  faster  decay  of  the 
excess-vibrational  energy  with  smaller  system  size.  The  results  in  frames  (c)  and  (d) 
seem  to  support  this  fact  with  the  results  from  the  smaller  bath  simulations  relaxing  to  a 
much  lower  final  excess-vibrational  energy.  The  HO2  embedded  in  125  Ar  atoms  is  ~4% 
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above  the  equilibrated  limit,  see  Fig.  19,  while  the  HO2  embedded  in  the  1000  Ar  is 
~10%  above  this  limit.  We  also  fit  the  curves  in  Fig.  22d  to  Eq.  46  and  the  fitting  results 
are  tabulated  in  Table  XL 
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Figure  22.  Relaxation  results  as  a  function  of  increasing  numbers  of  Ar  bath  atoms,  (a) 
Scaled  bath  gas  temperature;  250  Ar,  purple  curve;  500  Ar,  green  curve;  750  Ar,  blue 
curve;  1000  Ar,  yellow  curve,  (b)  Same  as  (a),  except  scaled  bath  pressure,  (c) 
Vibrational  energy  decay  on  semi-log  scale  with  the  same  coloration  as  in  (a),  dashed, 
black  lines  correspond  to  fits  to  a  single  exponential  function,  (d)  Vibrational  decay  with 
corresponding  fits  to  Eq.  46. 
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Table  XI.  Parameters  and  error  bars  as  a  function  of  number  of  Ar  atoms. 


#  of  Ar 

A 

B 

oab 

ki 

0  kl 

k2 

0  k2 

125 

0.26 

0.74 

0.05 

0.0301 

0.0075 

0.0058 

0.0003 

250 

0.26 

0.74 

0.05 

0.0253 

0.0053 

0.0049 

0.0003 

500 

0.25 

0.75 

0.05 

0.0229 

0.0051 

0.0042 

0.0002 

750 

0.28 

0.72 

0.05 

0.0201 

0.0040 

0.0039 

0.0002 

1000 

0.31 

0.69 

0.05 

0.0203 

0.0035 

0.0036 

0.0002 

The  amplitudes  of  the  exponentials,  red  curves  in  Fig.  23a,  appear  to  be  only  weakly 
affected  by  the  number  of  bath  atoms  within  the  calculated  error  bars.  The  magnitudes  of 
the  time  constants,  blue  curves,  show  a  decrease  with  increasing  number  of  bath  gas 
atoms.  These  results  indicate  that  the  mechanism  of  the  relaxation  process  is  only  very 
weakly  dependent  on  the  system  size,  but  the  speed  of  the  process  is  accelerated  by  the 
heating  of  the  bath  gas,  even  after  the  estimated  equipartition  energy  of  HO2  is  subtracted 
off.  The  magnitudes  of  both  time  constants  are  plotted  in  Fig.  23b  on  an  expanded  scale. 
The  larger  time  constant  increases  by  about  30%  going  from  the  largest  number  of  bath 
atoms  to  the  smallest;  however,  it  should  be  noted  that  the  calculated  error  bars  indicate 
there  is  a  great  deal  of  overlap  among  the  larger  time  constants.  The  smaller  time 
constant  increases  by  ~40%  going  from  the  largest  to  the  smallest  number  of  bath  gas 
atoms.  There  is  not  as  much  overlap  of  the  error  bars  of  the  smaller  time  constants, 
except  with  the  larger  simulation  size  results.  This  indicates  that  the  slower  process  is 
more  sensitive  to  the  issue  of  the  finite  size  of  the  simulation. 
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Figure  23,  The  variation  of  the  time  eonstants  and  amplitudes  obtained  in  fitting  Eq.  46 
as  a  funetion  of  the  number  of  bath  gas  atoms  at  800  K  and  1300  atm.  The  lines  are 
intended  to  guide  the  eye.  The  results  from  the  125  Ar  atoms  at  1300  atm  have  also  been 
included  in  the  plotted  data,  (a)  The  solid-blue  line  correspond  to  the  magnitude  larger 
exponential  and  the  solid-red  line  corresponds  to  the  associated  amplitudes.  The 
dashed-blue  line  corresponds  to  the  magnitude  of  smaller  exponential  and  the  dashed-red 
lines  correspond  to  the  associated  amplitudes,  (b)  Plots  of  the  magnitude  of  the  time 
constants  on  an  expanded  scale;  larger  time  constants  are  in  blue  and  are  relative  to  the 
left  y-axis  and  the  smaller  time  constants  are  in  yellow  and  are  relative  to  the  right  axis 
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D.  Conclusions 

The  relaxation  of  HO2  embedded  in  a  dense  Ar  bath  was  presented  in  this  ehapter. 
Our  work  was  motivated  by  the  nonlinear,  density  dependent  effects  on  relaxation 
reported  in  the  experimental  studies  of  Schwarzer  et  on  azulene,  and  Benzler  et 
on  cycloheptatriene,  and  in  the  simulations  done  by  Paul  et  on 

hexaflurobenzene.  Given  the  relative  importance  of  HO2  in  combustion,  we  were 
interested  in  seeing  whether  the  relaxation  of  vibrationally  excited  HO2  has  a  density 
(pressure)  dependence.  We  chose  to  simulate  the  relaxation  process  of  excited  HO2 
embedded  in  125  Ar  atoms  at  800K  and  35,  70,  100,  500,  900  and  1300  atm.  We  have 
established  and  validated  a  general  scheme  for  these  simulations  that  is  similar  to  the 
methods  proposed  by  others. 

We  found  that  the  decay  of  excess  vibrational  energy  was  well  represented  by  a  bi¬ 
exponential  function.  This  finding  is  surprisingly  similar  to  the  results  of  Paul  et  al, 
who  reported  simulations  of  a  much  larger  solute  in  a  diatomic  bath  gas.  We  also 
observed  that  both  our  time  constants  and  the  time  constants  reported  by  of  Paul  et  al. 
have  a  non-linear  densities  dependence  similar  to  the  behavior  reported  by  Schwarzer  et 
al.  and  Benzler  et  al.  We  also  reanalyzed  the  results  of  Paul  et  al.  and  found  similarities 
with  our  results  and  evidence  of  turnover  behavior.  We  found  that  our  lower  pressure 
results  do  not  extrapolate  to  zero  in  the  limit  of  an  infinitely  dilute  gas. 

We  also  investigated  the  effect  of  finite  sizes  on  the  relaxation  process.  The  primary 
size  effect  is  heating  of  the  bath  due  to  energy  lost  by  the  solute.  We  found  as  expected, 
that  the  bath  heats  up  to  temperatures  predicted  by  energy  equipartition.  This  heating 
does  not  appear  to  change  the  relaxation  mechanism,  but  rather  it  only  appears  to 
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accelerate  it.  The  effeet  of  heating  is  not  completely  accounted  for  by  making  the 
asymptotic  thermal  energy  of  HO2  a  function  of  the  temperature  of  the  bath  gas.  The 
finding  that  the  low-pressure  time  eonstants  do  not  extrapolate  to  zero  in  the  limit  of  low 
density  eould  be  due  to  this  finite  size  effect,  and  the  lower  pressure  time  constants  would 
likely  deerease  with  inereasing  system  size.  The  number  of  Ar  atoms  will  have  to  be 
increased  by  a  factor  of  eight  to  recover  the  low-pressure  limit. 

While  we  have  found  evidence  of  bi-exponential  relaxation  timescales  and  pressure 
effeets  on  the  relaxation  time  constants  we  do  not  have  a  elear  meehanism  for  what  is 
oeeurring.  It  is  probable  that  there  is  an  interplay  of  multiple  factors,  which  ultimately 
gives  rise  to  the  results  we  report  in  this  paper.  A  probable  meehanism  that  could  have 
two  relaxation  timescales  is  the  following.  We  hypothesize  that  the  first  timeseale  is 
related  to  the  relaxation  of  rotational  energy.  Large  extensions  of  the  0-0  bond  will 
eause  large  modulation  of  the  rotational  energy.  The  oscillation  of  the  0-0  bond  is 
intimately  tied  to  the  loss  or  gain  of  rotational  energy.  A  shorter  0-0  bond  length  means 
smaller  moments  of  inertia,  so  rotational  energy  will  be  more  easily  ehanged  by  collisions 
with  the  bath.  Longer  extensions  of  the  0-0  bond  would  mean  larger  moments  of  inertia 
that  would  make  the  radieal  are  more  resistant  to  changes  in  the  rotational  kinetic 
energy..  It  is  probable  that  loss  of  rotational  energy  is  mitigated  with  a  slow  leakage  of 
energy  from  the  0-0  vibrational  degree  of  freedom.  A  easeade  of  energy  into  the  bath 
via  the  rotational  degrees  of  freedom  would  eontinue  until  these  degrees  of  freedom  are 
thermalized  and  the  magnitude  of  vibrational  energy  is  not  sufficient  to  sustain  longer 
extension  of  the  0-0  bond.  A  potential  test  would  be  to  measure  the  range  of  oseillation 
of  0-0  as  the  radieal  is  progressively  loses  energy  by  eollisions.  The  relaxation  proeess 
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would  then  transition  to  the  slower  meehanism  which  would  likely  consist  of  a  slow 
direct  energy  transfer  out  of  the  lowest  frequency  vibrational  degrees  of  freedom. 

The  turnover  of  the  relaxation  time  constants  with  increasing  density  implies  a  change 
in  how  the  bath  gas  atoms  are  interacting  with  the  solute.  We  suspect  a  change  in  the 
local  environment  of  the  solute  could  be  the  cause  of  this  turnover  in  the  relaxation  time 
constants.  In  the  next  section,  we  discuss  some  future  directions  of  research  that  could 
helpful  in  understanding  the  bi-exponential  timescales  of  relaxation  and  how  the  bath  gas 
density  could  cause  the  turnover  of  the  relaxation  time  constants. 

E.  Directions  of  Future  Research 

This  initial  work  was  motivated  by  evidence  of  a  breakdown  in  the  IBC  assumption 
of  linear  scaling  of  relaxation  time  constants  with  increased  pressure.  While  our 
simulation  results  show  a  breakdown  in  the  expected  linear  scaling  effects  of  pressure  the 
mechanism  of  this  breakdown  has  not  been  resolved  and  there  are  several  avenues  for 
future  investigation.  It  has  been  suggested^^’’’*^^  that  a  clustering  effect  of  bath  gas  around 
the  solute  might  inhibit  collisional  energy  transfer  and  could  be  the  source  of  the  turnover 
effect.  Collision  events  can  have  a  varying  finite  lifetime,  so  it  seems  probable  that 
collision  events  could  overlap  at  higher  densities.  Ohmine  proposed  a  potential 
energy  criterion  for  determining  the  beginning  and  ending  of  collision  events,  and  a 
relative  scale  for  classifying  the  collisions  that  were  most  effective  in  relaxing  an  excited 
solute.  Collisions  were  determined  to  occur  when  the  sum  of  the  intermolecular 
interactions  of  the  bath  with  the  solute  became  greater  than  some  threshold.  This  scheme 
only  identifies  collisions  from  the  cumulative  interatomic  potential  energy  without 
distinguishing  the  number  bath  gas  atoms  that  could  be  participating  in  a  collision  event. 
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A  more  precise  definition  would  look  at  changes  in  the  magnitude  of  individual 
interactions.  The  beginning  and  end  of  a  collision  could  then  be  determined,  which 
would  allow  assignment  of  the  duration  of  a  collision.  The  distributions  of  the  lifetimes 
of  the  solute -bath  collision  complexes  and  the  average  number  of  local  collision  partners 
could  clarify  how  many-body  effects  or  clustering  influence  density-dependent 
relaxation. 

The  decay  of  excess-vibrational  energy  was  well  represented  by  a  bi-exponential 
function  across  all  of  the  simulated  pressures.  As  previously  discussed  in  Sec.  C.2a,  a 
similar  finding  was  reported  by  Paul  et  for  a  hexaflurobenzene.  The  bi-exponential 
fitting  is  evidence  of  two  relaxation  processes.  It  is  possible  that  IVR  is  the  cause  of  this 
binary  timescale  behavior.  We  hypothesize  that  the  low  frequency  mode  (0-0)  is  more 
likely  to  lose  energy  by  collisions  than  high  frequency  modes  (0-H  stretch  and  H-0-0 
bending  motion).  The  energy  flow  out  of  the  0-0  mode  could  be  mediated  by  V-R 
coupling,  with  the  thermalization  of  the  rotational  degrees  of  freedom  occurring  relatively 
quickly.  The  high  frequency  modes  would  likely  remain  relatively  isolated  from  the 
effects  of  collisions  with  the  bath  gas.  A  slow  transfer  of  energy  out  of  the  high 
frequency  modes  could  repopulate  the  0-0  mode  causing  the  collisional  energy  loss  to 
repeat  until  the  chemical  species  is  thermalized.  Once  sufficient  vibrational  energy  is 
lost,  the  repopulation  of  energy  into  the  0-0  mode  would  slow;  consequently,  the  energy 
transfer  to  the  bath  gas  would  reflect  this  change  in  mechanism.  Our  IVR  simulations  of 
the  isolated  radical  show  the  possibility  of  slow,  picosecond  timescales  for  internal 
energy  transfer,  which  would  fit  in  this  mechanistic  scheme.  We  are  not  aware  of  any 
internal  energy  transfer  results  for  hexafiuorobenzene.  However,  von  Benten  et  al. 
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reported  pieoseeond  IVR  timescales  in  the  experimental  relaxation  studies  of  mono- 
substituted  benzenes,  so  it  does  not  seem  unreasonable  to  propose  a  slow  IVR  process  as 
the  source  of  the  bi-exponential  relaxation  reported  by  Paul  et  al. 
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V.  Conclusions 

Some  studies  of  the  dynamics  of  the  hydrogen  peroxyl  radical  have  been  reported  in 
the  previous  chapters.  Beginning  with  the  dynamics  of  the  isolated  radical  we  have 
compared  three  PESs  for  HO2  in  terms  of  phase  space  structures  using 
surfaces-of-section,  and  for  IVR,  isomerization,  and  dissociation  trajectory  dynamics. 
The  IVR  shows  multiple  timescales  with  the  relaxation  of  excited  00  stretch  being 
especially  slow  for  both  the  IMLS  and  XXZLG  PESs.  There  is  evidence  that 
isomerization  helps  facilitate  IVR  from  the  low  frequency  0-0  mode  to  the  higher 
frequency  modes  involving  the  motion  of  the  hydrogen  atom.  Slow  IVR  in  turn  leads  to 
prominent  mode-specific  effects  in  isomerization  (on  both  PESs)  and  dissociation  (only 
accessible  on  XXZEG).  The  isomerization  and  dissociation  populations  are  best 
described  as  bi-exponential  with  the  two  timescales  being  attributed  to  slow  IVR.  The 
DMBE  IV  results  contrast  with  the  results  obtained  with  the  other  surfaces:  IVR  occurs 
on  much  shorter  timescales,  mode  specific  effects  are  minor  or  nonexistent,  and  any 
effect  of  isomerization  on  IVR  is  also  minor.  The  unimolecular  reactions  involving  this 
surface  are  well  described  by  a  single  timescale. 

Euture  directions  of  study  on  isolated  HO2  could  be  devoted  to  further  quantification 
of  the  non-statistical  reaction  dynamics  by  comparison  of  results  from  additional 
trajectory  studies  and  variational  transition  state  theory.  Investigation  of  the  interaction 
between  rotation  and  vibration  would  also  be  useful  in  further  understanding  the 
timescales  of  energy  transfer  in  this  radical.  Comparative  studies  of  the  dynamics  of  the 
HO2  radical  and  the  hydrogen  peroxyl  anion  could  offer  additional  insights  into  the 
factors  and  timescales  of  the  dynamics. 


130 


147 


The  second  study  focused  on  the  relaxation  of  HO2  in  a  dense  gas  environment.  The 
relaxation  of  the  radical  by  collisions  with  an  Ar  bath  was  found  to  be  bi-exponential. 
The  relaxation  time  constants  were  found  to  have  a  nonlinear  density  dependence.  Finite 
size  effects,  primarily  in  the  form  of  heating,  do  not  appear  to  change  the  mechanism,  but 
rather  speed  up  the  rate  of  thermalization.  Such  finite  effects  are  a  likely  source  for  the 
imperfect  match  of  the  time  constants  to  the  limit  of  zero  density.  The  issue  of  finite  size 
would  likely  have  the  same  effects  on  the  convergence  of  the  lower  pressure  relaxation 
time  constants. 

Future  research  directions  are  proposed  with  emphasis  on  understanding  if  the 
suggested  clustering  effect  or  many  body  collisions  with  the  bath  could  be  the  source  of 
the  nonlinear  density  dependence  of  the  relaxation  results.  Ohmine  proposed  a  method 
for  identifying  and  quantifying  collisions,  which  we  have  modified  for  analyzing 
individual  collision  events.  We  suggest  that  IVR  and  clustering  of  the  bath  around  the 
solute  as  possible  sources  for  the  bi-exponential  decay  of  the  excess-vibrational  energy. 
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APPENDICES 
Appendix  I 

The  contents  of  this  appendix  is  a  listing  of  the  modifications  and  additions  that  were 
made  to  the  classical  dynamics  code  GenDyn  in  the  course  of  this  research.  The  core 
program  and  user  manual  were  obtained  from; 
http://www.chem.missouri.edu/Thompson/research/gendyn.htm 

The  subroutines  are  loosely  grouped  according  to  the  role  of  the  subroutine  in  the 
program:  reading  runtime  conditions,  selecting  initial  conditions,  integration,  potential 
energy  surface,  analysis.  The  user  manual  of  the  original  code  is  sufficiently  detailed  for 
a  straightforward  compilation  and  execution  of  the  code.  The  modularity  of  the 
programming  style  allows  for  easy  inclusion  of  new  routines  and  algorithms.  The 
majority  of  the  following  routines  were  programmed  during  the  course  of  this 
dissertation,  unless  indicated  with  italics.  The  subroutines  in  italics  are  generally 
composed  of  code  that  was  primarily  written  by  others.  These  subroutines  were 
interfaced  and  modified  for  use  within  GenDyn  and  original  authorship  is  acknowledged 
in  line  with  the  code.  Some  of  the  modified  routines  are  core  routines  from  the  original 
GenDyn  code  where  multiple  authors  likely  contributed;  these  codes  are  not  in  italics. 
Much  of  the  usage  of  the  code  is  described  in  Ch.  Ill  and  IV  and  liberal  comments  in  line 
with  the  code  are  present  throughout  the  subroutines. 

Some  of  the  code  functionality  of  the  routines  related  to  the  dense  gas  simulations 
were  not  used,  but  could  be  of  future  interest.  The  additional  functionalities  that  were  not 
used  were:  periodic  boundary  conditions;  interatomic  force  fields  for  Lennard-Jones 
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interactions  in  addition  to  Buckingham  potential;  and  interactions  of  a  solute  embedded 
in  either  diatomic,  or  triatomic  bath  gases. 

Read  in:  cellinput.f,  cellinput.pbc.f,  eqcartmod.f,  inputmod.f,  varymod.f, 

Initial  conditions:  canongas.f,  cellboltz.f,  cellboltz.pbc.f,  celladj.f,  cellsz.f,  cellsz.pbc.f, 
diag.f,  hessgen.f,  inmnrg.f,  randno.f,  stmom.f,  roto.excit.f, paxis3.f90 

Integration:  diffverlet.f,  diffvprod.f,  diffadm.f,  pbc.f 

PES:  dmbeiv.f,  dynamic _parameters.f90,  genpotgas.f,  genpotgas.pbc.f,  genpotmod.f,  , 
ho2gxpot.f,  ho2mb2001f  ho2pesfast.f,  intem.pbc.f,  utility./,  eqnmotmod.f, 
PE S _3D  subroutine. fdO, 

Analysis:  analysis. f,  backinteg.f,  eckart.f,  hamsys.f,  press. f,  rotation.f,  rdf  f,  writeout.f 
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Appendix  II 

Supplemental  Material  for  “A  Classical  Trajectory  Study  of  the  Intramolecular  Dynamics 
and  Unimolecular  Dissociation  of  HO2”  by  J.  Perry,  R.  Dawes,  A.  F.  Wagner,  and  D.  L. 
Thompson,  J.  Chem.  Phys.  139,  084319  (2013). 

This  supplemental  material  concerns  (1)  details  of  the  XXZLG  potential,  (2)  kinetic 
model  for  energy  transfer,  and  (3)  the  determination  of  latency  time. 

A.  XXZLG  details 

In  the  XXZLG  routine  any  input  geometry  is  converted  to  Rqo,  Roh,  and  ZOOH  and 
then  splines  are  fit  to  the  ab  initio  energies  in  these  coordinates.  Our  initial  trajectories 
and  subsequent  analysis  of  the  XXZLG  PES  revealed  two  minor  issues  that  can  be 
loosely  characterized  as  inflexible  coordinates  and  a  deficient  database. 

The  coordinate  inflexibility  is  a  consequence  of  how  XXZLG  forces  exchange 
symmetry  with  a  valence  bond  coordinate  system  defined  above  that  does  not  display  that 
symmetry.  The  shortest  0-H  bond  is  always  selected  for  spline  interpolation.  This 
choice  has  the  desired  effect  that  any  path  through  the  0-0  perpendicular  bisector  plane 
will  be  continuous  and  exhibit  the  required  exchange  symmetry;  however,  exchange 
symmetry  also  requires  that  any  perpendicular  path  through  the  bisector  plane  have  zero 
slope  at  the  intersection  of  the  plane.  The  XXZLG  routine  does  not  enforce  this 
requirement,  resulting  in  a  discontinuous  slope  at  the  perpendicular  bisector  plane.  The 
patch  discussed  in  Ch.  Ill  remedies  this  problem. 

The  XXZLG  database  for  the  spline  fitting  is  deficient  in  two  ways.  First,  it  would 
appear  that  at  least  for  collinear  geometries,  many  of  the  electronic  structure  calculations 


134 


151 


are  not  fully  converged  to  the  ground  state.  Within  the  database  there  are  173  pairs  of 
collinear  geometries  that  correspond  to  the  same  Cartesian  geometry;  e.g., 
(Roh,^OOH,Roo)  =  (1.0ao,180°,2.2ao)  or  (3.2ao,0°,2.2ao).  The  two  energies  for  each  pair 
should  be  identical  but  the  root-mean-square  (rms)  difference  between  the  two  energies 
over  the  173  pairs  is  ~3.3  kcal/mol.  For  energies  below  70  kcahmol,  which  encompasses 
both  the  H  +  O2  and  O  +  OH  asymptotes,  there  are  22  pairs  with  a  1.4  kcal/mol  rms 
difference  and  a  maximum  difference  of  4.2  kcal/mol. 

A  second  deficiency  in  the  database  is  its  sparsity  of  configurations.  There  are  two 
measures  of  this.  The  first  can  be  illustrated  for  the  isomerization  saddle  point.  The 
closest  database  geometry  to  that  saddle  point  is  (Roh,^OOH,Roo)  =  (2.2ao,50°,2.7ao). 
In  evaluating  the  energy  for  (2.2ao,ZOOH,2.7ao),  XXZLG  will  evaluate  the  spline  shown 
in  Fig.  1  where  the  dots  are  the  available  database  entries  for  ZOOH  given  (Roh,Roo)  = 
(2.2ao,2.7ao). 
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Fig,  1  Spline  of  the  XXZLG  potential  energy  with  respeet  to  ZOOH  for  database  entries 
where  (Roh,Roo)  is  fixed  at  (2.2ao,2.7ao).  The  inset  is  the  region  of  the  saddle  point. 

Those  entries  are  for  10°  angular  inerements  and  feature  -300  keal/mol  energy 
ehanges  over  a  single  inerement  whieh  leads  to  ~1  keal/mol  ripples  in  the  spline  near  40° 
and  50°  as  seen  in  the  insert  of  Fig.  1.  Using  the  eleetronie  strueture  methods  that  were 
used  to  eompute  the  XXZLG  data,  we  ealeulated  an  energy  1 .3  keal/mol  higher  than  the 
XXZLG  energy  at  the  XXZLG  isomerization  saddle  point  geometry,  eonsistent  with  the 
ripple  features  in  Fig.  1.  A  seeond  measure  of  database  sparsity  arises  from  the  faet  that 
~25%  of  the  geometries  in  the  database  have  Ron  values  that  are  the  larger  of  the  two 
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Roh  values  for  the  given  (Roo,^OOH).  Because  XXZLG  only  uses  the  ab  initio 
database  at  the  smaller  of  the  two  possible  Rqh  values,  XXZLG  will  not  exactly  recover 
the  database  energy  for  these  geometries.  About  25%  of  the  database  is  not  fully  used  or 
put  another  way,  25%  of  the  database  can  be  used  to  quantify  the  XXZLG  “fitting  error.” 
For  energies  less  than  70  kcal/mol,  the  rms  fitting  error  is  0.17  kcahmol  for  ~1,300 
geometries.  The  convergence  errors  for  collinear  geometries  discussed  above  contribute 
to  this  rms  error,  but  when  these  points  are  removed  the  rms  error  is  0.13  kcal/mol  with  a 
maximum  error  of  1 .0  kcal/mol. 

Of  course  any  PES  based  on  an  ab  initio  method  will  have  errors  inherent  in  the 
specific  electronic  structure  theory.  However,  the  XXZLG  PES  has  fitting  errors  that  are 
measured  in  a  few  tenths  of  a  kcal/mol  globally  and  a  few  kcal/mol  for  selected  regions 
for  energies  under  70  kcal/mol.  These  fitting  errors  can  only  be  improved  by  more 
complete  convergence  of  the  electronic  structure  theory  calculations,  especially  in  the 
collinear  region.  Eurther  improvements  would  require  either  a  denser  database  or  a 
different  fitting  method,  such  as  interpolative  moving  least  squares,^^‘^  that  can  be 
supported  by  a  database  of  non-uniform  point  distribution  constructed  to  minimize  the 
fitting  error. 
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B,  Kinetics  Model 


As  a  means  for  interpreting  the  energy  transfer  in  isolated  HO2  Albert  F.  Wagner 
proposed  and  we  jointly  developed  an  analytical  kinetics  model  of  energy  transfer  in  a 
three  mode  system  that  is  based  on  two  principles:  (1)  the  change  in  the  energy  of  a  mode 
is  directly  proportional  to  the  amount  of  energy  the  mode  currently  has  and  (2)  the  energy 
of  each  mode  asymptotically  relaxes  to  the  same  value,  i.e.,  the  total  energy  divided  by 
the  number  of  modes.  These  two  principles  produce  the  following  kinetics  model: 
dE, 


dt 

dE^ 

dt 

dE^ 

dt 


~~  (^12  ^13)-^!  ^12^2  ^13-^3 


—  ki2Ei  {k^2  +  ^23)-^2  ^13-^3 


^  ^12-^1  ^12-^2  (^13  ^23)-^3 


(S.l) 


where  Ej  is  the  energy  in  mode  i  and  kij  is  the  IVR  rate  constant  for  the  transfer  of  energy 
from  mode  i  to  mode  j.  The  second  principle,  (2)  above,  of  asymptotic  equipartition  of 
energy  requires  kij  =  kp.  Equation  (1)  does  not  make  explicit  the  conservation  of  total 
energy.  By  a  change  of  variables  from  {Ej,  E2,  E3)  to  {E1+E2+E3  =  E,  Ei  -  E2  =  A,  E2  - 
E3=  E),  becomes  the  fixed  total  energy  while  A  and  /"asymptotically  go  to  zero  to 
satisfy  the  equipartition  requirement.  Applying  the  change  of  variables  to  Eq.  (1)  results 
in: 


dt 

dA 

—  =  -(2^12  +  ^13)^  -  (^13  -  ^23)r 

dt 

dT 

— —  =  —(2^23  +  kjjjr  —  (kjj  —  ky2)A 

dt 


(S.2) 
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The  analytic  solution  for  is  =  Etotai  and  the  solutions  for  A  and  F  are 
A(t)  = 

.  r  3  ^  z.  1  (S.3) 

r(0  =  -cAe-'^'‘  +  —  — A{t) 

W)  2 

where  ki(2)  =  ki2+k]3+k23  ±  A,  c  =  {4A  -  2  kn  -  ku  -  k23)l{ki3-k23),  A  =  (ku^ +ki3^ +k23^  - 

1  /9 

knkn  -  ki2k23  -  kj3k23)  and  A  and  B  are  as  yet  undetermined  constants.  Equation  (3) 
has  three  nonlinear  parameters  that  together  give  the  three  kij.  From  the  above,  Eqs.  (S.3) 
can  be  rearranged  to  provide  explicit  expressions  for  E/t)  : 

fl  fl  1  A 

‘  -r-P  U  3;  [2  4c  J 

fl  fl  1  A 

E,  (t)  =  E  - -  Ae-^'‘ - +  —  Be-^^‘  (S.4) 

v2  yz  4cJ 

(2r\  (  1  ^ 

EM)  =  E  +  —  Ae^’^'" - Be’^^' 

asymp  ^  2c  j 

Here  Easymp  is  the  asymptotic  final  energy  equal  to  one  third  of  the  total  energy. 

This  model  can  display  two  different  kinds  of  structural  deficiencies.  First,  as  results 
in  CH.  Ill  C.l  show,  optimized  values  of  ki,  k2,  and  c  can  be  reasonable  while  an 
underlying  ky  might  have  a  negative  value.  Finder  such  circumstances,  the  kinetics  model 
may  be  phenomenologically  useful  but  its  ky’s  are  not  capable  of  general  physical 
interpretation.  In  these  cases,  the  reported  ki,  k2,  and  c  values  are  for  a  constrained 
optimization  where  any  optimized  negative  ky  is  constrained  to  be  zero.  One  way  this 
behavior  can  happen  is  when  the  model  is  fit  to  IVR  decay  over  too  short  a  period  of  time 
relative  to  the  decay  constant  such  that  the  decay  is  very  nearly  linear  in  time.  Since  the 
energies  in  all  three  modes  depend  only  on  A(f)  and  r(t),  under  such  circumstances  the 
data  can  only  strongly  support  four  constants,  namely  two  intercepts  and  two  slopes  5a 


139 


156 


and  Sr.  However,  the  model  has  five  adjustable  parameters:  ki,  k2,  c.  A,  and  B  [see  Eq. 
(S.4)].  If  A  and  B  represent  the  intercept  information,  the  ki,  k2,  and  c,  or  more 
specifically,  three  ki/s  must  be  correlated  to  obtain  the  two  optimum  slopes.  This  can 
lead  to  a  ky  being  negative  because  large  excursions  including  into  negative  values  of  any 
one  ky  result  in  negligible  change  in  fitting  error  if  the  two  other  ky's  change  in  a 
correlated  fashion.  The  two  slopes  sa  and  sr  contain  all  the  rate  information  and  in 
CH.  Ill  C.I  some  of  the  discussion  will  involve  their  value. 

A  second  deficiency  arises  from  the  nature  of  Eq.  (S.l)  where  each  E,  is  controlled  by 
one  of  the  three  pairs  of  rate  constants  that  can  be  formed  from  the  three  %’s.  In  the 
extreme  case  where  two  rates  constants  are  zero,  then  one  of  the  E/’s  and  E  itself  are  time 
independent,  allowing  the  kinetics  model  to  display  two  different  asymptotes.  CH.  Ill 
C.I  does  have  such  a  case. 

The  parameters  A  and  B  can  be  made  to  give  exactly  the  initial  conditions  or  they  can, 
along  with  the  nonlinear  parameters,  be  least-squares  adjusted  to  the  trajectory  results  for 
A  and  E.  As  explained  in  CH.  Ill  C.I,  H  and  B  are  best  determined  by  a  linear  least- 
squares  process,  while  the  %  parameters  are  determined  by  standard  nonlinear  least 
squares  fitting  methods.  Altogether,  this  makes  five  adjustable  parameters  to  represent 
a  particular  IVR  simulation.  These  nonlinear  methods  were  complemented  by  a 
discretized  search  for  apparent  local  minima  on  various  3D  grids  in  %  parameter  space. 
This  discretization  approach  revealed  topological  features  of  the  rms  error  surface  with 
respect  to  the  %  parameters  that  were  useful  in  determining  the  global  minimum. 
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C,  Latency  time 

For  some  mode-specific  excitations,  for  an  initial  period  of  time  there  is  either  no 
decay  rate  or  an  accelerating  decay  rate.  When  decay  curves  display  such  features, 
population  data  over  an  initial  period  of  time  were  excluded  from  the  fit  of  Eq.  (20).  This 
exclusion  or,  in  effect,  latency  time  to  was  determined  by  contrasting  the  rms  fitting  error 
to  the  rms  sampling  error.  The  rms  sampling  error  comes  from  the  fact  that  only  a 
finite  number  of  trajectories  are  used  to  construct  the  population,  implying  that  the 
population  at  each  time  has  a  sampling  uncertainty  whose  square  summed  over  every 
time  increment  leads  to  the  rms  sampling  error.  Both  the  fitting  error  and  the  sampling 
error  change  as  more  data  is  excluded  from  their  calculation  by  an  increasing  latency 
time.  However,  the  fitting  error  is  initially  larger  than  the  sampling  error  and  rapidly 
decreases  as  latency  time  increases,  all  because  of  the  unsuitability  of  Eq.  (20)  for  fitting 
early  time  features.  At  the  same  time,  the  sampling  error  changes  little  with  latency  time 
because  excluded  sampling  errors  are  small  when  the  normalized  population  is  near 
unity,  which  it  is  at  early  times.  We  define  the  latency  time  as  the  earliest  time  at  which 
the  fitting  and  sampling  rms  errors  are  equal. 
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